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ABSTRACT 


A  method  of  analyzing  laminated-composite,  linear-elastic  plate  and 
shell  structures  has  been  developed  based  on  the  hybrid  stress  finite-element 
model.  Flat-plate  elements  are  used  for  both  plate  and  shell  structures. 

Three  quadrilateral  elements  designated  by  ELEMl,  ELEMZ ,  and  ELEMR,  all  with 
linear  in-plane  boundary  displacement  and  linear  in-plane  stress  assumptions, 
are  chosen  for  use.  Transverse  shear  deformation  effects  and  rotary  as  well 
as  in-plane  inertia  are  included  in  all  of  these  elements.  The  ELEMl  element 
can  be  applied  to  single-layer  thin  or  moderately  thick  plates  and  shells. 

The  ELEMR  element  can  be  applied  to  multilayer  thin  plates  and  shells,  and  the 
ELEMZ  element  may  be  used  to  analyze  single-layer  or  multilayer,  thin  or  thick 
plates  and  shells. 

The  stiffness  matrix,  mass  matrix,  and  equivalent  thermal  loading  vector 
of  these  elements  that  are  used  in  static,  vibration,  and  thermal  stress 
analyses,  respectively,  are  developed  and  applied  to  analyze  various  plate  and 
shell  problems.  Comparisons  of  results  are  made,  whenever  possible,  with  other 
existing  solutions.  This  present  method  of  analysis  provides  accurate,  effi¬ 
cient  predictions  of  (a)  displacements  and  stresses  under  static  or  thermal 
loading  and  (b)  natural  frequencies  under  free  vibration.  The  method  developed 
in  this  study  can  be  ised  as  a  reliable  tool  in  structural  analysis  and  design; 
it  can  also  be  applied  to  explore  material  property  characterizations  and 
selections  for  laminated  composite  materials. 
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SECTION  1 
INTRODUCTION 


1.1  Purpose  of  the  Study 

♦ 

In  the  design  of  structures  consisting  of  layered-fiber  reinforced  com¬ 
posite  materials,  a  designer  is  often  faced,  with  the  problem  of  a  complex, 
structure  subjected  to  complex  loading  conditions.  The  structure  is  usually 
in  the  form  of  combined  shells  and/or  plates.  The  loadings  may  include  static, 
vibratory,  and/or  transient  mechanical  loading  afe  well  as  thermal  loads,  i  In 
order  to  assure  that  the  structure  can  sustain  various  loading  conditions  > 
without  inducing  failure  of  its  normal  functions,  the  designer  must  choose  the 
most  appropriate  composite  material,  and  determine  the  appropriate  dimensions 
and  configuration  of  the  layered  composite  material.  Therefore,  reliable 
analytical  tools  are  needed  (1)  to  predict  the  behavior  of  a  tenti  live  design 
and  (2)  to  improve,  if  necessary,  a  itentative  design  based  on  detailed  analysis 
More  specifically,  the  designer  should  be  able  to  predict  (1)  the  stress  and/or 
strain  distribution  at  any  point  within  the  layered!  composite  and  (2)  deflec¬ 
tions  of  the  structure,  so  that  the  failure  or  safety  of  the  structure  can  be 
predicted  (a)  by  comparing  the  stress  or  strain  against  an  appropriate  strength 
failure  criterion  and/or  (b)  by  comparing  the  calculated  deflection  against 
specified  deflection  limits.  If  the  tentative  design  is  not  satisfactory , 
then  the  analysis  will  point  cut  the  weak  points  of  the  structure  and  suggest 
changes  in  the  design,  which  could  mean  both  a  change  in  dimension  and  in 
orientation  of  the  fiber  directions  of  the  composite. 

The  main  objective  of  this  project  is  to  develop  appropriate  structural 
analysis  tools  to  assist  the  designer  in  material  selection  and  material  charac 
terization.  In  particular,  it  is  intended  that  these  analysis  tools  when  ade¬ 
quately  developed  and  validated  be  employed  to  make  a  realistic  assessment  of 
the  types  and  values  of  composite  mater  ials/structurer.  properties  that  will 
provide  "optimum"  resistance  to  the  transient  mechanical  and/or  thermal  loads 
*  experienced  by  the  structure.  Because  o-f  the  complex  nature  of  the  problem, 
the  finite-element  method  is  selected  as  an  appropriate  and  versatile 
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mathematical  model  to  cope  with  the  complex  structure  and  loading  conditions. 
Because  transverse  shear  deformation  effects  are  known  to  be  important  in 
laminated  composites  11,2,3,4)*,  the  usual  classical  (Kirchhoff)  thin  plate 
and  shell  theory  cannot  be  used  directly;  accordingly,  the  finite-element 
.model  developed  and  utilized  in  the  present  study  includes  transverse  shear 
deformation  effects.  Since  the  conventional  displacement  finite- element 
iqethod,  when  applied  to  relatively  thick  laminated  plates,  either  has  failed 
to  predict  accurately  the  local  deformations  and  stresses  of  a  plate  under 
bending  [5]  or  is  too  expensive  to  use  because  too  many  degrees-of- freedom 
are  involved  for  even  relatively  simple  problems  (6) ,  the  hybrid  stress 
finite- element  model  pioneered  by  Pian  (7)  provides  a  promising  alternative. 

,  The  hybrid  stress  model  of  the  finite-element  analysis  is  characterized 
by  an  assumed  stress  field  in  the  interior  of  the  element  as  well  as  an 
; assumed  displacement  field  along  the  boundary  of  the  element.  Since  a  stress 
field  is  directly  assumed,  the  prediction  of  stresses  in  a  given  problem  is 
generally  more  accurate  than  that  of  the  conventional  assumed-displacement 
model  where  only  the  displacement  field  is  assumed  for  both  the  interior  and 
along  the  boundary  of  the  element.  Also  great  flexibility  in  terms  of  the 
formulation  and  the  selection  of  element  features  is  achieved  because  the 
assumed  interior  stress  distribution  and  the  assumed  boundary  displacement 
distribution  functional  can  be  selected  independent  of  each  other.  The  com¬ 
patibility  of  the  interelement  displacements  is  easily  satisfied,  since  only 
the  boundary  displacement  is  assumed,  while  in  the  conventional -displacement 
model  the  satisfaction  of  compatibility  usually  poses  a  difficult  problem. 

More  important  is  the  fact  that  transverse  shear  deformation  effects,  which 

cannot  be  neglected  in  many  laminated  composite  plate  and  shell  problems,  can 

1 

be  easily  taken  into  account  by  the  hybrid  stress  model  without  an  increase  in 
the  number  of  element  degrees-of- freedom,  as  indicated  in  Ref.  S3.  Because  of 
transverse  shear  deformation  effects,  thick  laminated  plates  may  have  different 
cross-sectional  rotations  for  each  layer.  Thus,  a  severe  warping  of  the 
original  plane  section  of  the  whole  laminate  develops  under  certain  loading 

«  ~ 

References  are  denoted  by  numbers  in  square  [  ]  brackets. 
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conditions.  This  phenomenon  can  be  closely  modeled  by  the  hybrid  stress  model 
with  ease  as  will  be  seen  later  in  this  report. 

Structures  of  prime  interest  for  future  analysis  and  design  include, 
among  others,  aerospace  vehicles  in  the  shape  of  a  low-aspect-ratio  lifting 
oody,  irregular  shaped  shells,  and/or  axisymmetric  shapes.  Cutouts,  branches, 
stiffeners,  etc.  may  be  present.  Layered  anisotropic  construction  may  comprise 
a  significant  portion  of  the  structure. 

For  the  analysis  of  curved  shells,  much  recent  finite-element  work  has 
been  carried  out  using  curved  finite  triangular  and/or  quadrilateral  elements 
for  cylindrical  shells  [9-11],  conical  shells  [12],  and  shells  of  revolution 
[13,14].  These  finite  elements  are  convenient  if  the  structure  is  smooth  and 
regular  since  they  can  model  the  actual  structure  faithfully  geometrically; 
however,  these  elements  are  not  suitable  or  convenient  for  analyzing  shells 
of  irregular  shape  (irregular  boundaries,  arbitrarily-oriented  cutouts, 
arbitrarily  curved  snape;  ...  )  Also,  it  should  be  noted  that  curved-shell 
elements  are  developed  for  some  specific  restricted  geometry  and  hence  when 
used  for  modeling  irregular  curved  shells  would  entail  significant  geometric 
modeling  errors  unless  the  element  size  were  very  small  compared  with  the  main 
dimensions  of  the  shell.  Further,  if  the  element  size  must  necessarily  be 
very  small,  it  would  be  adequate  to  model  the  curved-shell  structure  by 
flat-plate  elements.  In  this  case  the  analyst  enjoys  the  advantages  of  (a) 
simplicity  in  element  stiffness  derivation  and  computation  and  (b)  great 
computer  program  versatility  in  being  able  to  treat  a  wide  variety  of  irregu¬ 
larly-shaped  structures  with  relatively  simple  computer  program  logic  and  com¬ 
putations  since  the  element  compatibility  conditions  and  the  transformations 
from  local  element  coordinates  to  a  global  coordinate  system  for  the  entire 
structure  are  handled  in  a  uniform  fashion. 

Flat-plate  elements  based  upon  the  assumed  displacement  approach  have 
been  used  to  analyze  curved  shells  [15].  Recently,  Wolf  [10]  used  flat-plate 
hybrid  stress  finite  elements  to  analyze  the  "pinched  cylindrical  shell";  his 
predictions  compare  favorably  with  an  independent  solution  obtained  with 
assumed  displacement-type  curved  finite  elements  [13].  Further,  the  assumed 
stress  hybrid  model  nas  been  shown  to  be  applicable  readily  to  plates  which 
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involve  transverse  shear  deformation  and/or  thermal  loading  [dj . 

based  upon  the  above  considerations  and  especially  its  relevance  to  the 
present  objectives  of  this  project,  the  hybrid-stress  finite-element  model  has 
been  selecteu  for  use  in  this  study.  In  particular,  flat-plate  elements  are 
developed  for  both  plate  and  shell  problems. 

Before  any  problems  involving  complex  transient  loading  conditions  can 
be  solved,  it  is  essential  to  develop  and  establish  the  reliability  of  the 
pertinent  finite  elements  under  static  loading  conditions,  bike  conventional 
displacement  finite-element  models,  the  hybrid  stress  model  for  static  analysis 
is  also  characterized  by  a  stiffness  matrix.  Once  the  stiffness  matrix  of  an 
element  is  proved  to  be  suitable  for  use,  then  dynamic  problems  can  be  pursued 
by  developing  a  mass  matrix  and  by  first  applying  it  to  vibration  analyses  to 
establish  the  reliability  of  tne  mass  marrix. 

As  for  thermal  problems,  the  most  important  task  is  to  develop  an  equiva¬ 
lent  loading  vector  corresponding  to  a  given  temperature  distribution. 

Therefore,  in  the  past  twelve-month  period,  effort  has  been  concentrated 
on  the  selection  and  verification  of  the  behavior  of  selected  finite  elements 
under  static  loading  conditions,  on  the  development  and  verification  of  the 
mass  matrix  of  each  of  several  finite  elements,  and  on  an  equivalent  thermal 
loading  vector. 

1.2  Synopsis  of  the  Investigation 

A  method  of  analyzing  plate  and  shell  structures  built  with  layered- 
fiber-reinforced  composite  materials  has  been  developed.  The  work  done  in  the 
past  twelve  months  is  described  in  the  following  sections  nf  this  report. 

In  Section  2,  the  static  analysis  aspect  is  described.  It  begins  with 
a  derivation  of  the  stiffness  matrix  for  single-layer  flat-plate  elements. 

Then,  several  candidate  elements  are  studied  and  compared  against  each  other 
on  two  test  examples.  It  is  found  that  triangular-based  quadrilateral  elements 
are  inferior  to  "basic"  quadrilateral  elements;  therefore,  a  basic  type  of 
quadrilateral  element  is  selected  as  the  "best"  element  to  be  used  in  this 
study.  This  single-layer  element  is  extended  to  a  multilayer  element  in 
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Subsection  2.3  which  also  includes  a  general  derivation  of  the  stiffness  matrix 
for  a  multilayer  element.  Both  single-layer  and  multilayer  elements  are  tested 
on  plate  problems  in  Subsection  2.4.  Example  problems  include  single-layer 
isotropic  plate  problems.  Both  thin  plates  and  thick  plates  are  included.  In 
Subsection  2.5,  the  flat-plate  elements  are  applied  to  shell  problems,  in¬ 
cluding  shells  of  revolution.  The  necessary  transformation  for  srngle-layor 
anu  multilayer  flat-plate  elements,  when  applied  to  shell  geometries,  are 
developed  in  Subsections  2.5.1  and  2.5.2,  respectively.  Example  problems  are 
discussed  in  Subsection  2.6,  where  problems  of  thin  cylindrical  shells  and 
conical  shells  of  different  thicknesses  are  treated.  In  both  the  plate  and 
shell  examples,  comparisons  are  made  with  exact  elasticity  solutions  or  other 
approximate  solutions  whenever  possible.  Excellent  accuracy  is  observed.  Pre¬ 
dicted  stresses  are  sliown  to  be  within  2%  of  the  exact  elasticity  solution  in 
one  case.  These  examples  are  discussed  in  Subsection  2.7. 

Section  3  begins  with  a  discussion  on  the  variational  basis  for  selecting 
a  mass  matrix  in  the  hybrid  stress  model.  Then  a  hybrid-rational  mass  matrix  as 
well  as  a  lumped  mass  matrix  is  developed  and  compared.  It  is  found  that  the 
hybrid-rational  mass  matrix  is  slightly  better  than  the  lumped  mass  matrix  for 
thin  plate  vibrations.  Since  it  is  expected  that  the  former  will  be  much  better 
than  the  latter  for  thicker  plate  vibration  analyses,  the  hybrid-rational  mass 
matrix  is  used  thereafter  throughout  this  study.  Example  problems  of  plate  and 
shell  vibration  are  given  at  the  end  of  Section  3.  The  predictions  for  the 
lowest  and  next  higher  mode  frequencies  are  compared  with  other  solution  method 
predictions  and  are  found  to  be  very  accurate. 

Thermal-stress  analysis  is  described  in  Section  4.  For  given  tempera¬ 
ture  distributions,  an  approximation  of  the  temperature  within  an  element  is 
first  established.  Then  an  equivalent  loading  vector  is  calculated  based  on 
this  approximate  representation  of  the  actual  temperature  distribution.  Since 
numerical  examples  of  thermal  stress  analysis  do  not  appear  frequently  in  the 
literature,  in  most  of  the  examples  given  in  this  section,  comparisons  cannot 
be  made.  However,  the  equivalent  loading  vector  is  verified  through  simple 
examples. 


A  summary  of  the  work  for  this  period  of  study  and  concluding  remarks 
are  yiven  in  Section  5.  This  report  ends  with  an  appendix  on  the  programming 
aspect  of  the  present  finite-element  analysis. 


6 


SECTION  2 
STATIC  ANALYSIS 

2.1 Formulation  of  Single-Layer  Flat-Plate  Elements 

2.1.1  Deri  vation  of  Element  S  tif  fj  a  os  s  Matrices 

Tne  formulation  of  the  element  stiffness  matrix  by  the  hybrid  stress 
model  is  based  on  the  Principle  of  Minimum  Complementary  Energy  {Kef.  7) . 

Along  each  boundary  of  the  element,  an  assumed  misplacement  distribution  is 
selected;  these  element-boundary  displacements  are  then  expressed  in  terms  of 
the  nodal  displacements  q.  hext,  an  assumed  stress  distribution  in  terras  of 
undetermined  stress  parameters  b  is  chosen  tnroughout  the  interior  of  the  ele¬ 
ment;  this  same  distribution  function,  if  uesired,  may  also  be  used  to  approxi¬ 
mate  tne  prescribed  surface  tractions.  In  terms  of  these  distribution  func¬ 
tions  and  associated  unknowns  q  and  £,  one  may  write  an  expression  for  the 

complementary  energy  if  .  By  setting  dir  /3b  =  0  to  minimize  the  complementary 

c  c  ~ 

energy,  one  obtains  an  approximate  solution  for  b_  in  terms  of  q.  This  enables 
one  to  evaluate  the  internal  strain  energy  U  entirely  in  terms  of  the  gen¬ 
eralized  nodal  displacement  unknowns  q.  Then,  by  recognizing  tiiat  the  stiff- 

T 

ness  matrix  appears  in  the  following  form  for  U  as  U  =  1/2  q  k_  q,  one  there¬ 
by  identifies  k  in  the  hybrid  stress  model  formulation.  This  procedure  is 
presented  in  detail  in  the  following. 

The  total  complementary  energy  of  one  finite  element  is  given  by 

TTc  =  ~2  j v  -S  ijitfi  <Tij  CTtf  d  1/  T;  U;  6  k  (2.1) 

where 

S.  a  compliance  constants  of  the  material 
ljki 

0 .  .  =  stresses 

id 

V  =  volume 

A  =  boundary  area 

T.  =  surface  traction 

l 

u  =  prescribed  boundary  displacement 
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The  functional  n  is  a  minimum  when  the  stresses  satisfy  the  equilibrium  con- 
c 

dition. 

For  convenience,  consider  a  statically  loaded  continuum  which  is  repre¬ 
sented  by  a  single  finite  element.  In  deriving  the  element  stiffness  matrix, 
the  displacements  along  each  boundary  of  the  finite  element  are  expressed  in 
terms  of  the  nodal  displacements  q  and  certain  interpolation  functions  L, 
such  that  the  displacement  compatibility  conditions  with  the  neighbor ing  ele¬ 
ments  (when  one  models  the  continuum  with  a  number  of  finite  elements)  are 


satisfied , 


U  = 


(2.2) 


The  element  stresses  in  the  interior  of  the  element  are  then  expanded  in  terms 
of-  a  finite  number  of  stress  parameters  £ 

sr=e£  (2.11 

where  P  is  chosen  to  satisfy  the  homogeneous  equilibrium  equations.  The 
surface  tractions  can  then  be  written  in  the  form  of 

I  =  £B  (2.4) 

where  R  is  obtained  by  applying  the  element  boundary  conditions  to  P. 
Substituting  Eqs.  2.2  to  2.4  in  Eq.  2.1,  one  obtains, 

^  -  z  €  <1  f s  e  ■ ^  -  f  (i  sTi  m  )  t 


By  defining 


H  =  j  PTSPdv 

§  =S,  STtdA 


(2.6) 


Eq.  2.S  becomes 


rrc  =  {t 


The  best  approximate  solution  for  8  for  the  problem  is  obtained  by  setting 


dTT  /jp  to  zero.  Tlie  result  is 
e  ~ 
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H  £  “  £  g  -  o 

from  which, 

p  =  y'G  g 


(2.y) 


The  first  term  in  Eq.  2.7  for  7i  is  the  total  strain  energy  1)  in  the  element. 


By  substituting  Eq.  2.9  into  U  one  obtains 

U  =  j  STaViS  5- 

/C  »w  ■*-  /%»  /w 


By  definition 


U  =  T  % 


(2.10) 


(2.11) 


0  0  ^  v  '  •  ■L-L' 

-O  -v 

where  is  the  element  stiffness  matrix. 

Thus,  by  the  hybrid  stress  model,  the  element  stiffness  matrix  is  given 

by 

is  =  si  y  ^  (2.i2) 

2.1.2  Material,  Stress,  and  uispl acement^ _As_s umptio n s 

The  matrix  ]5  in  Eq.  2.6  is  the  compliance  matrix  that  relates  stresses 
£  to  strains 


€  =  5  <T 

^  it-  — 


(2.13) 


In  this  st.uay,  the  following  square  symmetric  material  properties  are 
assumed  for  each  layer  of  the  composite  material: 


.  r\  \ 

/ 

^  _ .. 

&rr 

e%H 

I  \  1  »  I  - 

»  _  ru  _  " 

Ei  Ez  E3 

iht  I  _  jjjj 

El  E 2  E3 

± 'jl  l>31  I 

E,  Ez  £3 


1  \ 
|0x  I 


(2.14a) 


1  p-— 


i*^.ih!AM**p 


m 


with  E2  =  E3,  v12  =  V13.  V21  =  v31.  V23  =  vy2,  U13  =  G12,  and  "1"  being  the 
longitudinal  fiber  direction. 

In  case  the  material  axes  do  not  coincide  witli  the  element  coordinate 
axes  (see  Fig.  la) ,  the  following  transformed  £  is  used: 


'$„m?*sun4 

+(iS,i+S„)tnxn' 


5«(rrt<+’iV+ 

(5„ 

Saire+U 


o 


*(ZSli+Ss')rf7)1 

\ 

\ 

V Mum1  su 

\ 

O  °  Surt+Sa*1 


o 

(2$rlSa'Ssr)*h 


x 

0  0 

( ~1  $ii  ~$u)  2^]*4r^  ^ 

*Wl 


x  4(S,-lS/i  +-%.-25rj.)wV 
^  +  (w\  jft’Srf 

> 


(2.1<lb) 


where  m  =  cos  a,  n  =  sin  a,  anu  =  the  corresponding  constants  in  Eq.  2.14b. 

The  in-plane  interior  (and  boundary)  stress  field  P  for  0  ,0  ,  and  0 
- c -  ~l  x  y  xy 

consists  of  a  stretching  part  P  and  a  bending  part  P  : 


P,  =  Ps  (*,5)  +  Pj> 


(2,15a) 


In  this  study,  P  is  assumed  to  be  linear  in  x  and  y,  and  PL  is  either  linear 
~s  ~b 

in  x  and  y  or  quadratic  in  x  and  y.  The  out-of-plane  stress  field  £2  for 

0,0,  and  0  is  determined  from  the  equilibrium  equations  of  elasticity, 
xz  yz  z 

The  matrix  £  in  Eq.  2.3  consists  of  and  P^.  Different  assumptions  for 
different  elements  are  summarized  in  Table  1  and  a  detailed  description  can  be 
found  m  Subsection  2.2. 
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Along  element  boundaries  (Fig.  11/),  the  following  displacement  field 
is  assumed,  for  a  typical  side  1-2  of  an  element: 

U  =  U,<l-S)  +  UiS  +  Zfeyi(l-S)f9^Sj 

v  =  v,  ( i-s)  +  uz 5  -  z  f  e>,  r/'5;t  sj 

=  W,  (l-S)  +•  1^2  5  (2.16) 

where  s  is  the  coordinate  measured  along  the  side  1-2,  and  nodal  displacements 
u,  v,  w,  0^,  and  0^  are  defined  in  Fig.  lb.  In  other  words,  the  in-plane  dis¬ 
placements  u  and  v  are  linear  in  both  z  and  s,  while  the  transverse  displace¬ 
ment  w  is  linear  in  s  but  is  independent  of  z.  The  number  of  degrees  of  free¬ 
dom  per  node  is  5.  Both  quadrilateral*  and  triangular  elements  are  considered 
in  this  study  (Figs,  lc  and  Id);  therefore,  the  total  degrees  of  freedom  per 
element  is,  respectively,  either  20  or  15. 

Variations  of  Eq.  2.16  can  be  obtained  by  adding  a  mid-point  node  and 
changing  the  linear  functions  of  s  in  Eq.  2.16  to  quadratic  functions. 

Based  on  the  assumed  interior  (and  boundary)  stress  field  Pb  in  Eq.  2.15 
and  the  assumed  element  boundary  displacement  field  Lq  in  Eq.  2.16,  the  element 
stiffness  matrix  k  can  be  obtained  by  carrying  out  the  appropriate  integration 
of  Eq.  2.6  and  substituting  into  Eq.  2.12. 

It  should  be  noted  that  shear  deformation  effects  are  taken  into  account 

in  Eq.  2.16,  where  0  and  0  are  the  rotations  of  the  normal  element  of  a  plate; 

x  y 

these  quantities  are  not  to  be  identified  with  the  slopes  w  and  w  of  the 

» y 

middle  surface.  In  other  words,  the  normal  to  the  undeforroed  middle  surface 
dues  not  remain  normai  to  tne  deformed  middle  surface.  Also,  in  the  evaluation 
of  H  by  Eq.  2.6,  the  transverse  shear  energy  terms  are  included. 

Various  elements  with  different  assumed  P3  and  Lq  are  studied  and  the 
results  are  described  in  the  following  subsection. 

2. 2  Parametric  Study  and  Selection  of  Single-Layer^  Elements 

Since  the  stress  field  of  an  element  can  be  assumed  independently  of  the 
boundary  displacement  field,  it  is  possible  to  match  different  pairs  of 

* 

Included  are  square,  rectangular,  and  irregular  quadrilateral  geometries. 
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assumptions  and  to  select  those  elements  with  the  better  convergence  behavior. 

The  stresses  in  an  individual  lamina  may  be  considered  as  the  combina¬ 
tion  of  the  stretching  and  bending  components  as  shown  in  Eq.  2.15.  When  the 
stretching  component  is  represented  by  a  complete  linear  expansion  in  x  and  y, 
while  the  bending  component  is  represented  by  a  complete  quadratic  expansion 
in  x  and  y,  the  stress  field  can  be  represented  by  34  independent  stress 
parameters  as  follows: 

!  2  3  ♦  S  l  r  8  9  fo  il  12  13  h  IS  Ifc  17  18  I?  20  2/  K  ZJ  2*  25  2(,  27  ?8  2?  3o  J I  ft  33  H 


In  constructing  the  various  models,  one  of  the  following  three  differ¬ 
ent  stress  fields  is  used  in  a  giver  model: 

(a)  Type  (L,L  ):  This  corresponds  to  columns  1-‘J,  10,  11,  12,  15, 

10,  iy,  and  20  of  the  matrix  in  Eq.  2.17,  with  a  total  of  10  B's. 
The  designation  (L,L  )  stands  for  a  complete  linear  expansion  in 
x  and  y  for  the  stretching  component  and  an  incomplete  linear  ex¬ 
pansion  m  x  and  y  for  the  bending  component. 

(b)  Type  (L, L) :  This  corresponds  to  columns  1-20,  with  a  total  of 
20  p's-  The  designation  (L,L)  stands  for  a  complete  linear 
expansion  in  x  and  y  for  both  the  stretching  and  the  beniing 
component. 


(c)  Type  (L,y)  •.  This  corresponds  to  columns  1-34,  with  a  total  of 

34  p's.  The  designation  (L,y)  stands  Cor  a  complete  linear  ex¬ 
pansion  in  the  stretching  components  and  a  complete  quaciratic 
expansion  in  bending. 

Three  types  of  element-boundary  displacement  fields  (along  side  1-2, 
Cor  example)  are  used: 


(a)  Type  (L) :  This  corresponds  to  tne  following: 


u  =  u,  o-sj  +  UiS  +  zf^o-sj  +  e^si 
v  =  V,  U-S)+  v2  s  -Z[Qx,0-s)+Bz2s} 


(2. Id) 


W  =  IV,  (I-  5)+  W2S 

where  s  is  the  coordinate  along  one  side  of  the  element.  The 
designation  (L)  stands  for  a  complete  linear  expansion  in  s. 

(b)  Type  (Q) :  This  corresponds  to  the  following: 


U  -  u,  Ci-s)  +  UzS  +  4uma-s)S  +z[9yto-sj+6j2 s] 

ir  =  V,  u-s)  +  y2s-t4vmfi-s)5-z[0,r((/-s;t0;ri5]  {2  i9) 

WlU-5)^WzS*4Wm(i-S)$ 

The  subscript  in  stands  for  the  mid-point  side  node.  The  designa¬ 
tion  (Q)  stands  for  a  complete  quadratic  expansion  for  the  dis¬ 
placements. 

(c)  Type  (Q  ) :  This  corresponds  to  the  following: 
ll  —  «,(i-s;  +  «2S+Z(^,U-S)+ 

V  -  V,  (p-s)+l65-Z[  Qxld-s)+  9x*sl  (2.20) 

w  =  v,  a-s3+vJs  +  (enJ-©nj(f-s)si/? 
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wnere  0  is  the  rotation  along  the  element  side  and  is  a  function  of  0  and  0  . 

n  x  y 

The  designation  (Q  )  stands  for  a  quadratic  expansion  tor  w  only.  If  type  (Q) 
is  used  on  certain  sides  and  type  (L)  is  used  on  the  remainder,  then  the 
designation  is  (L,Q) . 

Because  shear  deformation  effects  are  of  major  concern  in  this  study, 
the  selected  elements  should  be  able  to  converge  rapidly  both  for  thin  and 
thick  plates.  Therefore,  two  extreme  test  problems  are  used  to  evaluate  '  ari- 
ous  elements  in  order  to  select  the  most  suitable  elements  for  future  use.  One 
is  a  simply- supported  square  thin  plate  with  a  thickness- to-side  ratio  of  1:100, 
and  the  other  is  a  simply-supported  square  thick  plate  with  a  thickness-to- 
sidc  ratio  of  1:4.  The  former  is  centrally  loaded  and  the  latter  is  uniformly 
loaded . 


A  brief  description  of  the  elements  studied  and  the  associated  example- 
problem  solutions  are  summarized  in  Table  1  and  in  Figs.  2  and  3.  Mot  every 
result  is  plotted  in  Figs.  2  and  3  because  some  are  far  below  or  above  the 
range  of  the  plots.  A  total  of  nine  elements  has  been  tried;  these  are 
described  in  the  following,  together  with  brief  assessment-type  comments: 


(a)  ELKtU:  General  quadrilateral  element  with  stress  type  ( L,L) ,  displace- 
ment  type  (L) . 


This  element  is  suitable  for  problems  of  plates  or  shells  of  revolution. 
Convergence  is  fast  and  computation  of  the  k  matrix  is  inexpensive.  The 
total  oegrees-of- freedom  per  element  is  20. 


(b)  ELEMA:  Triangular  element  with  stress  type  (L,L)_  and  displacement  type  (L) 
This  element  does  not  converge  to  correct  answers  within  reasonable  mesh 
refinement  limits.  This  element  is  too  "stiff"  and  is  rejected. 

(c)  ELEMB :  Triangular  element  with  stress  type  (L.L)  and  displacement  tyjje 

(b,Q) 

To  improve  ELEtIA,  a  midpoint  node  is  added  to  the  interior  boundary  (see 
Table  1) .  The  improvement  is  not  enough.  This  element  is  still  too 
stiff  and  is  rejected. 
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(d)  E1.EMC:  Triangular  element  with  stress  type  (L,_E)  and  displacement  type  (y) 

To  make  ELEMfl  more  flexible,  midpoint 'nodes  are  added  to  every  boundary. 
Unfortunately,  the  stiffness  matrix  becomes  singular  even  after  the  5 

usual  boundary  constraining.  This  is  due  to  the  fact  that  the  number  of 
unwanted  displacement  modes  whicl^  induce  no  stresses  within  the  element 
(these  inodes  are  called  kinematic  inoues)  are  too  many  (4  for  bending 

only)  to  be  completely  suppressed  by  the  imposed  boundary  conditions  > 

(Ref.  17).  This  element  is  rejected. 

(e)  EEEMU:  Triangular  element  with  stress_type  (L, L_J  and  displacement  type  (L,^)) 
Another  way  of  making  an  element  more  flexible  is  to  reduce  the  number  of 
b's.  From  EEEMB,  the  number  of  £'s  is  reduced  to  10.  This  element  con¬ 
verges  rapidly  as  can  be  seen  in  Figs.  2  and  2.  however,  its  application 

!  i 

to  future  thermal  analysis  is  limited  because  its  stress  field  is  too  sim¬ 
ple.  Therefore,  it  is  rejected.  i 

(f)  ELEliE:  Triangular  element  with  stress  type  (L,Q)  and  displacement  type  (L) 

To  suppress  the  unwanted  modes  of  ELEMC,  a'  complete  quadratic  bending 

stress  fielu  is  used.  Due  to;  the  increase  in  the  number  of  stress 

\  : 

parameters  ii,  i. t  is  too  stiff  and  is  rejected. 

i 

(y)  KLEl'll' :  Triangular  element  with  stress  type  (L,Q)  ^nd  displacement  type  (g) 
Midpoint  nodes  are  added  to  LEE ME  to  make  it  more  flexible.  It  turns  out 
that  EEEMF  is  good  for  thin  plates  but  is  too  flexible  for  ti^ick  plates 
as  can  be  seen  from  Figs.  2  and  3. 

()')  EEEMG :  Triangular  element  with  stress  type  (L,g)  and  displacement  type  (g.  ) 

i 

This  element  has  the  advantage  of  increasing  order  of  w  without  increas¬ 
ing  the  number  of  degrees-of- freedom  per  node.1  The  convergence  is  moderate. 
Since  only  w  is  quadratic,  this  element  is  not  suitable  for  ihell  problems 
where  u,  v,  and  w  should  be  consistent.  Therefore,  it  is  rejected. 

(i)  E EEL  12 ; _ Four  triangular  elements  wi_th  stress  type  (E,g)  and  displacement 

type  (L,g) 

Four  triangular  elements  assembled  into  one  quadrilateral  element.  The 
convergence  for  a  thick  plate. is  very  rapid  aiui  the  convergence  for  a 
thin  plate  is  fairly  rapid.  This  element  can  be  used  for,  general  shell 
problems . 

1L> 
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Note  that  for  solution,  one  quarter  of  the  plate  was  modeled  by  uniform-sized 
finite  elements.  The  number  of  degrees  of  freedom,  DOF,  indicated  in  Figs.  2 
and  3  denotes  the  DOF's  for  the  quarter  plate. 

Out  of  the  nine  elements  tested  (see  Table  1) ,  only  two  of  them  show 
fast  convergence.  The  triangular-based  quadrilateral  element,  ELEM2,  is  more 
sophisticated  than  the  quadrilateral  ELEMl.  Consequently,  more  computing  time 
is  needed  to  calculate  the  stiffness  matrix  of  ELEM2.  Since  ELEM2  shows  no 
improvement  over  ELEMl,  ELEMl  is  used  hereinafter  throughout  this  study.  In 
the  following,  ELEMl  is  extended  to  represent  a  multilayer  element. 


2.3  Formulation  for  Multilayer  Flat-Plate  E lements 

2.3.1  Derivation  of  Element  Stiffness  Matrices 

In  the  present  formulation  of  stiffness  matrices  for  laminated  plate 
finite  elements,  the  self-equilibrating  stresses  for  each  lamina  are  initially 
assumed  independently.  The  stress  equilibrium  conditions  for  the  stresses  at 
the  interlamina  boundaries  are  then  introduced  as  a  finite  number  of  conditions 
of  constraint 


A(1  =  O  (2.21) 

Using  the  method  of  Lagrange  Multipliers,  the  complementary  energy  functional 
of  Eq.  2.7  is  modified  to 


(2.22) 


Again,  by  setting  3tt  */<)S=0  for  each  3,  one  obtains 

c 

HQ-Cjf  +  ATA=0 

from  which 

<t =  y"&  s  -  h"ata 


(2.23) 


(2.24) 


Substituting  into  Eq.  2.21,  results  in 

ANh6§-AHhAtx=c> 


(2.25) 
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or 


6  =  (ANVrcAtffiL; 

Substituting  X  from  Eq.  2.26  into  Eq.  2.24  and  then  substituting  the  resulting 
3  into  the  strain  energy  U,  one  obtains 

u  =  j  f(C3Y§-§W  <&>?£)"  Ati'V  $  12.27, 

Thus,  the  stiffness  matrix  is  given  by 

K  =  -6THVAT(AH''ATMH^T  (2.28) 


2.3.2  Stress  and  Displacement Assumptions 

Multilayer  laminated  thick  plates  often  exhibit  a  severe  warping  of 

the  cross  section  as  a  result  of  the  transverse  shear  deformation  effects  and 

the  discontinuous  material  properties  from  layer  to  layer  {Refs.  1,  2,  3,  4). 

To  model  this  warping  phenomenon  closely,  the  rotational  degrees-of-f reedom  0 

and  0y  should  be  assumed  to  be  different  for  each  layer.  Within  each  layer, 

0  and  0  are  still  assumed  to  be  a  constant  in  the  transverse  direction;  thus, 
x  y 

0^  and  0  can  be  represented  by  the  inplane  displacements  u  and  v  at  the  ipper 
and  lower  surface  of  that  layer  (Fig.  4).  The  transverse  displacement  w  is 
still  assumed  to  be  a  constant  in  the  transverse  direr*- ion  for  all  layers. 
Therefore,  the  total  uegrees-of-freedom  at  a  node  consists  of  one  w  and  (n+1) 
inplane  displacements  u  and  v,  with  n  being  the  number  of  layers.  The  number 
of  degrees-of-f reedom  per  node  is  then  2{n--  I)+1. 


For  a  typical  layer  i,  the  displacements  along  side  1-2  are: 


v  "  =J-(v,  ♦iskz-s)- +/2  -v,)V-sh 


Wf  7  =  7V,  0-S)+  W2S 


(2.29) 
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For  stress  distributions,  an  independent  set  of  P.jj^  (the  same  as  that 
used  for  ELEM1)  should  also  be  assumed  for  each  layer.  The  interface  equi¬ 
librium  is  achieved  by  demanding  that 

at  the  top  of  the  ith  layer 

—  /(t  /r  )  at  the  bottom  of  the  (i+l)th 

—  t  7.^,0 layer  (2.30) 

Equation  2.30  constitutes  a  constraint  on  the  stress  parameters  and  can  be 
arranged  in  tnc  form  of  Eq.  2.21.  Then  the  stiffness  matrix  of  a  multilayer 

element  can  be  computed  using  Eq.  2.20.  This  new  element  is  designed  as 

ELEMZ . 

For  relatively  thin  laminated  plates  and  shells,  the  warping  phenomenon 
is  not  severe  and  the  usual  assumption  that  a  plane  cross  section  remains 
plane  is  still  applicable.  The  corresponding  change  in  the  element  stiffness 
matrix  can  be  easily  achieved  by  demanding  that  the  nodal  displacements  u  and 
v  be  linear  in  the  transverse  direction  and  be  replaced  by, respectively ,0 

and  just  as  in  Eq.  2.1G.  Then  the  total  number  of  degrccs-of- freedom  per 

node  is  back  to  5  again  and  tiiis  reduced  multilayer  clement  will  be  designated 
as  ELK MR. 


2.4  Flat-Plate  Evaluation  Examples 

Several  plate  problems  nave  been  solved  and  a  comparison  of  results  is 
made  with  respect  to  exact  solutions  or  other  enproximatc  solutions.  In  all 
but  two  examples,  convergence  studies  arc  made  to  establish  the  reliability 
of  the  elements  EEEM1 ,  ELEMR,  and  ELE'IZ.  Problem  descriptions  and  the  finite- 
element  modelings  usr-H  are  given  in  Tables  2,  3,  and  4  for,  respectively,  single¬ 
layer  plates  and  three-layer  plates. 

2.4.1  Single-Layer  Isotropic  Hates 

(a)  Simply-Supported  Square  Thin  Plate  (Table 2,  Fig.  2) 

This  is  a  problem  used  in  the  selection  of  elements.  The  errors 
of  the  predicted  central  transverse  displacements  wc  for  elements 
1,  2,  A,  U,  and  F  arc  plotted  for  different  meshes  in  a  quarter  of 


tno  plate;  advantage  is  taken  ot  symmetry  so  that  only  one-quarter 
of  tne  plate  is  modeled  by  equal-sized  quadr ilateral  elements,  The 
percentage  errors  are  calculated  by  comparison  with  tne  exact  solu¬ 
tion  of  Kef.  In.  From  Fig.  2,  it  is  seen  that  tiie  error  diminishes 
very  rapiuly  ana  tnat  MLLHl  provides  somevnat  superior  results. 

( b )  Clamped  :>< juare  Thin  Plate  (Tal tie 2 ,  Fig.  5) 

To  see  whether  or  not  boundary  conditions  affect  the  convergence 
of  ELEMl  predictions,  a  clamped-plate  problem  is  solved.  A  plot 
similar  to  tnat  of  Fig.  2  is  presented  as  Fig.  i  and  shows  fast 
convergence.  Comparison  is  based  on  tne  exact  solution  of  Ref.  Id. 

( c )  Clamped  Rectangular  Thin  Plate  (Table  2,  I’ig.  G) 

3oth  a  square  element  and  a  rectangular  element  (1:2  aspect  ratio) 
for  ELL’Ml  are  used  to  see  the  effects  of  a  large  aspect  ratio  of 
the  element.  The  results  arc  compared  with  the  exact  solution  of 
Ref.  Id.  From  Fig.  6,  it  is  observed  that  in  both  cases  tiie  con¬ 
vergence  is  fast  although  the  pattern  of  convergence  is  reversed. 

(d)  Skewed  Simply-Supported  Thin  Fla to  (Table  2,  Tig.  7) 

Since  ELEM1  is  a  general  quadrilateral  element,  it  is  not  re¬ 
stricted  to  rectangular  plate  problems.  Various  solutions  {19,20,21] 
of  tiie  sxowed  plate  problem  identified  in  Table  2  were  collected  in 
Ref.  10  and  are  reprouuced  here  in  Fig.  7  together  witn  tiie  present 
ELL/li  solutions.  It  is  seen  that  tiie  present  simple  element  provides 
solution  accuracy  and  efficiency  comparable  to  other  higner-order 
elements  even  though  the  free  traction  condition  at  the  edge  is  not 
enforced  in  the  present  solutions  (tnis  enforcement  would  require 
a  special  element) . 

(e)  Simply-Supported  Square  Tnick  1‘latc  (Table  2,  Fig. _ J_> 

This  is  also  a  problem  used  in  tiie  selection  of  elements.  Com¬ 
parison  of  tne  predicted  central  displacement  is  made  with  tile 
exact  solution  of  Ref.  22.  It  is  seen  tiiat  111,1211  converges  even 
faster  in  tain  case.  This  is  expected  since  l.Li-Jf'.l  is  designed  to 
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cope  with  shear  deformation  and  in  this  problem  the  shear  de¬ 
formation  is  indeed  important  (about  23%  of  the  total  deforma¬ 
tion)  . 

The  results  of  the  above  examples  demonstrate  tne  reliability  of  KLEM1 
for  isotropic  plates.  Since  the  present  stuuy  is  mainly  intended  for  com¬ 
posite  plates,  the  above  results  can  only  be  regarded  as  .a  good  start,  and 
furtner  evaluations  have  been  carried  out  as  discussed  in  the  following. 

2.4.2  Two- Layer  Plates 

I'our  problems  of  two-layer  plates  are  solved.  The  ELEMR  element  is 
used  for  the  first  three,  since  they  involve  very  thin  plates  with  negligible 
warping.  For  the  last  one,  which  is  a  thick  plate  problem,  ELEHZ  is  used. 

(a)  Clamped  Cross-Ply  Rectangular  Thin  Plate  (Table  3,  Fig.  8) 

A  clamped-edge  rectangular  plate  consisting  of  two  layers  of 
fiber-reinforced  composite  materials  with  orientations  of  0°/y0° 
with  respect  to  the  positive  x-axis  is  uniformly  loaded.  The 
central  deflection  is  calculated  using  ELEMR  for  the  same  plate 
with  or  without  the  assumption  of  negligible  shear  deformation 
effects;  since  the  plate  is  indeed  very  thin  (with  a  side-to- 
thickness  ratio  of  1:100),  the  two  present  deflection  predictions 
differ  by  only  1%  and  are  both  very  close  to  but  slightly  larger 
than  the  CRT  (classical  laminated  plate  theory)  approximate  solu¬ 
tions  of  Ref.  23. 

(b)  Clamped  Cross-Ply  Square  Thin  Plate  (Table  3,  Fig,  9) 

As  indicated  in  Fig.  0,  this  problem  is  solved  twice  using  ELEMR. 

In  the  first  solution,  the  wuole  plate  is  oriented  in  the  usual 
way,  but  in  the  second  solution,  it  is  rotaLed  43°  so  that  the 
material  axes  are  no  longer  parallel  to  the  reference  coordinate 
axes.  The  purpose  is  to  nake  sure  that  the  transformation  of  the 
compliance  matrix  nas  been  correctly  programmed,  since  both  solu¬ 
tions  should  be  tne  same  if  the  transformation  is  correct.  Indeed, 
these  two  solutions  are  identical  and  converge  very  rapidly  in 


comparison  with  the  CPT  solution  of  Ref.  23. 

(c)  Clamped  Angle-Ply  Square  Thin  Plate  (Table  3,  Fig.  10) 

The  problem  considered  is  similar  to  the  previous  one  but  the 
fiber  orientations  are  +  43°.  Again,  the  two  solutions  of 
LLKMR  are  identical  and  converge  rapidly;  that  is,  to  within  3% 
of  the  CPT  solution  of  Ref.  23.  It  should  be  noted  that  the 
CPT  solution  tends  to  give  less  accurate  predictions  for 
angle-ply  problems  because  tne  distortion  of  the  plate  is  more 
severe  than  that  of  the  cross-ply  problems.  Hence,  the  5%  dif¬ 
ference  here  can  be  attributed  largely  to  the  error  of  the  CPT 
solution. 

(d)  Simply-Supported  Angle-Ply  Square  Thick  Plate  (Table  3,  Pig.  11) 

To  investigate  the  shear  deformation  effects,  a  thick-plate  prob¬ 
lem  is  solved  using  ELEMZ  for  different  orientations  of  the  two- 
layer  angle-ply.  The  present  solutions  differ  from  the  approximate 
solutions  of  Ref.  24.  Again,  the  difference  can  be  attributed  to 
the  error  of  the  approximate  solution  since  other  cases  tested  in 
Ref.  24  are  also  less  accurate  for  thick  laminated  plates.  Because 
of  the  shear  deformation  effects,  the  central  deflections  are 
approximately  doubled  compared  with  that  for  the  case  of  zero 
transverse  shear  deformation.  It  should  be  pointed  out  that 
studies  of  this  kind  will  assist  the  designer  in  the  selection 
of  the  appropriate  material  orientations  to  fulfill  certain  specific 

Ut'Si'jii  oLijkC  liVC3. 

2.4.3 _  Three-layer  Plates 

( a )  Infinite  Cross-Ply  Plate  with  Sinusoidal  Loading 
(Table  4,  rigs.  12a-12d) 

The  problem  of  an  infinitely  long  tnick  plate  loaded  sinusoidally 
in  the  short  direction  and  simply-supported  along  the  two  edges  is 
a  two-dimensional  plane-strain  problem  and  can  be  solved  exactly 
according  to  elasticity  theory  as  shown  in  Ref.  1.  Using  the 
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present  finite-element  method  and  FLEiiZ,  this  problem  is  solved  by 
putting  a  uniform  mesh  in  one  strip  in  the  short  direction  since 
the  solutions  are  independent  of  the  long-direction  coordinate. 

The  displacement  and  stress  predictions  are  very  close  to  those 
of  the  exact  solution  as  can  be  seen  in  Figs.  12a  through  I2u.  It 
is  observed  that  the  other  assumed  displacement  finite-element  so¬ 
lution  {Ref.  5} ,  which  also  takes  transverse  shear  deformation  ef¬ 
fects  into  account:  (1)  does  not  give  accurate  detailed  predictions 
but  only  average  values  of  cross-section  rotation  and  normal 
stresses,  and  (2)  is  very  close  to  the  CPT  solution,  which  docs  not 
take  shear  effects  into  account  at  all. 

(b)  Simply -Supported  Cross-Ply  Plate  (Tables  4  and  5,  Fig.  13) 

This  finite-dimension  problem  is  solvea  exactly  in  Ref.  2. 

Two  examples  with  different  aspect  ratio  and  thickness  ratio 
are  given  here  (Table  4,  cases  1)1  and  b2;  both  are  solved  using 
LLEMZ .  The  first  one  is  a  moderately  thick  rectangular  plate 
under  bending;  the  maximum  displacement  and  stress  components 
predicted  by  the  method  are  within  2%  of  the  exact  solution  [21 
as  can  be  seen  in  Table  5.  The  second  example  is  a  thick  square 
plate  under  bending;  the  stress  and  displacement  distributions 
at  crucial  sections  are  plotted  in  Fig.  13.  Also  plotted  are 
solutions  obtained  by  using  a  3-dimensional  assumed-displacement 
finite  element  (6J;  it  is  seen  that  the  present  method  gives  prac¬ 
tically  the  same  accuracy  with  many  fewer  degrr  cs-of  freedom 
(225  DOF  vs.  990  DOF). 

2.5  Application  to  Shell  Problems 

2.5.1  Transformation  for  Single-Layer  Elements 

In  the  application  of  flat-plate  elements  to  shell  problems,  it  is 
necessary  to  transform  the  clement  nodal  degrees-of-freeuom ,  which  are  in 
element  coordinates,  to  a  common  or  global  nodal  degree-of- freedom  system.  As 
discussed  in  Ref.  25,  the  transformation  can  be  expressed  in  the  following  form 
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CostXJ)  C os  (X  T)  Co$(X,Z) 

c&str.x)  Coscy.YJ  g>s cr.z) 


kH  =  CosiZ.x)  CosCZ.Y;  0>s CZ,Z) 


FLfMOJr 


Cos(x.x)  G>S<X,r)  Ck&.Z) 

Qn 

Cos(r.x)  Gsfr.Y)  0>s(y.2)  Q 


(2.31a) 


S  =  I-  £ 


(2.31b) 


where  the  "barred"  system  is  the  common  (global)  nodal  system. 


If  the  z-axis  of  the  common  nodal  coordinate  system  is  chosen  to  be  in 
the  direction  of  the  outer  normal  of  the  shell  surface  (Ref.  15) ,  then  the 
contributions  of  cos(x,z)  and  cos(y,z)  in  the  last  column  are  small  compared 
with  other  terms  and  can  be  neglected.  The  resulting  transformation  is  a  5x5 
transformation  and  the  nodal  degrees-of-freedom  consist  of  only  u,  v,  w,  0^, 
and  C  : 

y 


U  CosiX.X)  CoS  CX .  Z ) 

V  Co5(Y,  x;  CosCY.Y)  fofYZ) 

fc/U  £os(Z,  X)  GjfZ.Y)  CosiZ.l) 


t 

V  . 

ELEMENT 


CvsW.v  cos(x. y;  Bg 

Ces(Y.X)  Co'jY-Y)  By 


t40De  (2.32a) 


*  =  T s.r  t 


(2.32b) 


The  corresponding  element  stiffness  matrix  is 


K  =  TT  fc  T 


(2.33) 


Once  k  is  obtained  from  Eq.  2.33,  this  transformed  element  stiffness 
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matrix  can  be  assembled  into  a  system  stiffness  matrix  K,  as  usual,  and  the 
remaining  procedure  is  the  same  as  that  for  plate  problems. 

The  above  transformation  is  applicable  to  ELEMl  and  ELEMR. 

2.5.2  Transformation  for  Multilayer  Elements 

Since  the  nodal  displacements  of  a  multilayer  element  are  completely 

represented  by  the  interlayer  displacements  u,  v,  and  transverse  displacement  w, 

instead  of  u,  v,  w,  0  ,  and  0  ,  the  former  should  be  expressed  in  terms  of  the 

x  y 

latter  before  a  transformation  can  be  carried  out.  This  is  accomplished  by 

defining  an  equivalent  set  of  nodal  degrees-cf-freedom,  which  consists  of  the 

displacements  u,  v,  and  w  at  the  nodal  mid-thickness  station  of  an  element  and 

the  rotations  (0  ) .  and  (0  ) .  of  each  layer  of  the  element  at  a  nodal  station 
xr  y  l 

(see  Fig.  14) . 

The  relation  between  this  new  set  of  degrees-of-freedom  and  the  old  one 
can  be  put  in  matrix  form  as 


(2.34a) 


t  =  Tso  ?. 


(2.34b) 


For  a  typical  three-layer  case  as  show,  in  Fig.  17,  the  matrix  T  i: 

*wrr( ) 
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(2.35) 

Then  the  transformation  of  this  new  local  system  to  the  common  (global) 
system  is 


or 

S  =  I*  I 


(2. 3Ga) 

(2.3Gb) 
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where  the  submatrices  T  ,  and  T_  ,  are  the  submatrices  of  T_  _  of  Eq.  2.31. 

—•3x3  -"2x3  ~5x6 

The  overall  transformation  matrix  is  found  by  combining  Dq.  2.34  and 


Eq.  2.36: 

£  =  I*e  £  =  ~he  ? 

or 

5?  =  Tie  Te% 

The  transformed  element  stiffness  matrix  k  is  then 


(2.37) 


js  =  il%  a  Tti  (2*38j 

The  above  transformation  is  applicable  to  ELEMZ. 


2.6  Shell  Evaluation  Examples 

Thin  cylindrical  shell  problems  are  solved  first.  Then  conical  shel  s 
of  various  thickness  are  solved  to  demonstrate  the  effect  of  transverse  shear 
deformation.  The  problem  data  and  the  finite-element  modeling  employed  are 
given  in  Tables  6  and  7. 

2.6.1  Thin  Cylindr ical  Shells 

(a)  Pinched  Cylindrical  Shell  (Table  6,  case  (a) ;  Fig.  15) 

both  analytical  and  finite-element  solutions  [16,20,26]  are 
available  for  this  problem  and  are  collected  by  Wolf  in  Ref.  16. 
Although  the  present  solution  using  ELEMl  does  not  enforee  the 
free  traction  condition  on  the  simply-supported  edges,  its  ac¬ 
curacy  is  still  comparable  to  other  solutions  as  can  be  seen  in 
Fig.  15. 

( b )  Cylindrical  Shell  under  Ring  Load  (Table  6  c ass  (b) ;  Fig.  16) 

An  analytical  solution  for  this  problem  is  given  in  Ref.  27. 

Since  the  ring  load  is  axisymmetric ,  this  problem  can  be  re¬ 
duced  to  a  one-dimensional  problem  and  an  extensive  convergence 
study  can  be  performed  thriftily.  It  is  seen  in  Fig.  16  that  the 
deflection  under  the  ring  load  as  predicted  with  ELEMl  flat-plate 
elements  converges  rapidly  to  within  6%  of  that  of  the  analytical 
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solution.  An  increase  in  the  number  of  elements  eventually 
further  reduces  the  error  to  within  2t. 

2.6.2  Conical  Shells 

To  investigate  the  transverse  shear  deformation  effects  on  a  shell, 
several  conical  shell  problems  are  solved  with  ELEM1  and  ELEKR  flat-plate  ele¬ 
ments.  They  all  have  the  same  loading  and  geometry,  except  that  they  differ 
in  the  thickness  and  material  constituents.  The  chosen  loading  is  axisymmetric 
in  order  to  minimize  the  computing  cost. 

( a )  Isotropic  Conical  Shells  (Table  7,  Case  (a) ;  Figs .  17  to  20) 

First,  a  shell  with  thickness  =  0.025  in.  is  solved.  Com¬ 
parisons  of  displacement  and  moment  distributions  are  made 
in  Figs.  17a  and  17b  with  those  predicted  by  another  higher- 
order  finite-element  solution  {SABOR  4,  Ref.  28).  In  using 
45  elements  along  the  meridian,  the  present  solution  is,  within 
plotting  accuracy,  the  same  as  that  of  5ABOR  4.  The  latter  does 
not  take  transverse  shear  deformation  effects  into  account. 

Next,  the  thickness  is  increased  to  0.1  in.  (Figs. 18a  and  18b) 
and  then  to  0.5  in.  {Figs.  19  and  20).  There  are  still  no  ap¬ 
parent  transverse  shear  deformation  effects.  However,  when 
the  thickness  is  increase i  to  1  inch,  which  corresponds  to  a 
thickness-to- radius  ratio  =  1:15,  pronounced  transverse  shear 
deformation  effects,  reflected  by  the  difference  between  the 
present  solutions  of  ELEKl  and  those  of  the  SABOR  4,  are  ob¬ 
served  as  is  evident  from  Figs.  19  and  20-  For  a  thicker 
shell,  it  can  be  concluded  that  transverse  shear  deformation 
effects  must  be  taken  into  account  in  order  to  obtain  any 
meaningful  answer. 

(b)  Sandwich  Conical  Shell  (Table  7,  Case  (b)  ;  Fig.  21) 

Since  no  other  solutions  have  been  found  for  comparison,  only  a 
displacement  distribution  is  presented  for  this  example.  The 
sand v.i cli  sjiell  is  treated  as  a  three-layer  shell,  and  ELEMR  is 
used.  The  displacement  plot  in  Fig.  21  shows  a  rapid  die-out  of 
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the  loading  effect.  A  comparison  with  Fig.  17b  shows  that  the 
edge  deflection  of  the  sandwich  shell  is  smaller  than  that  of 
an  isotropic  shell,  which  uses  roughly  the  same  amount  of  materi¬ 
al  as  is  used  for  the  skin  of  the  sandwich  shell  (0.025  in. 
thick  vs.  0.02  in.  thick).  This  clearly  demonstrates  the  effec¬ 
tiveness  of  the  sandwich  structure  in  increasing  the  stiffness 
of  the  shell. 

( c )  Two-Layer  Conical  Shell  (Table  7,  Case  (c)  ■,  Figs.  20  and  22a-22b) 

To  simulate  some  missile-type  structures,  a  glass-phenolic  coating 
layer  is  placed  on  top  of  an  aluminum  base  layer  and  this  two-layer 
shell  is  analyzed  using  ELEMR.  The  normal  displacement  of  this 
two-layer  shell  is  plotted  in  Fig.  20  together  with  solutions  for 
other  similar  conical  shells.  Since  the  coating  layer  is  rela¬ 
tively  soft,  the  deflection  is  larger  than  that  of  a  similar  but 
isotropic  shell  of  1-in.  thickness.  The  maximum  meridional  normal 
stresses  in  both  layers  are  plotted  in  Fig.  22a.  As  expected,  the 
maximum  stresses  occur  at  the  two  outer  surfaces  of  the  shell. 

The  maximum  circumferential  normal  stresses  are  plotted  in  Fig.  22b. 
They  occur  at  the  interlayer  surface  and  at  the  inner  surface  of  the 
shell.  The  maximum  stress  in  the  base  material  is  more  than  three 
times  that  of  the  coating  material,  due  to  the  larger  stiffness  of 
the  base. 


It  is  clear  that  such  a  detailed  description  of  the  stress  distributions 
will  be  of  definite  value  to  a  designer  who  needs  this  basic  information  to  de- 
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2.7  Summary  and  Comments 


A  method  of  analyzing  fiber-reinforced  laminated  plate  and  she! 1  struc¬ 
tures  under  static  loading  has  been  developed.  The  hybrid  stress  finite-element 
method  has  been  chosen  as  the  mathematical  model  to  analyze  this  complex  problem, 
because  of  the  versatility  of  the  hybrid  stress  method,  problems  like  transverse 
shear  deformation  effects  and  multilayer  laminates  are  dealt  with  easily. 
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The  development  of  this  basic  analysis  tool-  has  been  carried  out  step 
by  step.  First,  several  candidate  sinyle-layer-plate  elements  have  been 
studied.  Based  on  two  test  problems,  a  quadrilateral  element  (ELEMl)  has  been 
chosen  for  later  use.  Then,  an  extension  to  multilayer  elements  (ELEMR  and 
ELEM2)  has  been  made,  and  various  plate  problems  have  been  solved  to  establish 
tiie  reliability  cf  the  present  method. 

Application  to  shell  problems  is  achieved , through  a  coordinate  trans¬ 
formation.  Then  several  shell  problems  have  been  solved  to  demonstrate  th6 
validity  of  this  transformation  and  to  study  transverse  shear  deformation  ef¬ 
fects  on  shells. 

In  addition  to  the  three  elements  selected  for  use  as  described  :in  this 
section,  another  similar  element  designated  as  ELEM3  was  developed  in  the  re¬ 
search  effort  under  this  contract  and  is  described  in  another  report  [29]*. 

That  element  is  designed  for  thin  multilayer  laminated  plates.  The  only  differ¬ 
ence  between  ELEM3  and  ELEMR  is  that  the  former  has  a  single  set  of  stress 
assumptions  for  tne  whole  flat  plate  element  reyardless  of  the  number  of  layers, 
while  the  latter  has  an  independent  stbess  assumption  for  eacn  layer  of  the 
flat-plate  element.  Studies  reported  in  Ref.  29  show  that  for  multilayer  plates 
with  tnickness-to-side  ratios  ranging  from  0.01  to  0.25,  ELEM3  provides  better 
results  than  does  ELFi.R.  The  errors  of  the  results  obtained  by  using  ELEM3  are 
uniformly  less  than  those  of  ELEMR.  however,  it  should  be  stressed  that  for 

.  i 

thick  laminates  with  a  typical  thiekness-to-side  ratio  equal  to  or  greater  than 
about  0.1,  the  multilayer  element  ELEM2  should  be  used  to  avoid  excessive 
error;  for  such  thickness  ratios,  ELE.M3  and  ELEMR  are  both  deficient. 


i 


In  Ref.  29,  botn  ELEMl  and  ELEMR  of  this  report  are  referred  to  as  ELEMl, 
since  they  nave  the  same  number  of  degrees-of-f reedom  per  element. 
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SECTION  3 


VIBRATION  ANALYSIS 

3.1  Vibration  Analysis  by  the  Hybrid  Stress  Model 

As  is  seen  in  the  previous  seetion,  the  static  hybrid  stress  finite- 
element  analysis  takes  the  same  form  as  does  the  displacement  finite-element 
analysis;  that  is,  each  element  is  represented  by  an  element  stiffness  matrix 
and  the  final  unknown  parameters  are  the  element  nodal  generalized  displace¬ 
ments.  In  other  words,  the  hybrid  stress  model  can  be  viewed  as  just  an  alter- 
•!  native  way  of  deriving  a  stiffness  matrix.  Along  this  line  of  thinking,  it  seens 
that  analogous  with,  but  somewhat  different  from  the  consistent  mass  matrix 
approach  used  in  the  displacement  finite  element  dynamic  analysis,  a  mass 
matrix  for  the  hybrid  stress  model  can  be  constructed  by  using  a  displacement 
field  for  the  interior  of  an  element  based  on  a  suitable  interpolation  of  the 

f 

assumed  boundary  displacements  of  the  hybrid  stress  element,  and  by  evaluating 
the  corresponding  kinetic  energy  of  the  element  in  terms  of  the  generalized 
nodal  velocities,  thereby  identifying  the  associated  mass  matrix.  This  was 
done  by  Dungar,  Severn,  and  Taylor  [30]  for  a  general  triangular  flat-plate 
element  and  a  right-angled  triangular  flat-plate  e’ement,  and  later  by  Dungar 
and  Severn  [31]  for  a  triangular  plate  element  with  variable  thickness.  Numeri¬ 
cal  examples  in  both  Ref.  30  and  Ref.  31  showed  convergence  of  the  predicted 
natural  frequencies.  This  approach  was  extended  to  treat  a  cylindrical  shell 
element  by  Hensnell,  Neale,  and  Warbuton  [32] ,  and  to  accommodate  a  triangular 
flat-plate  element  and  a  rectangular  flat-plate  element  by  Neale,  Henshell  and 
—  Edwards  133 J  - 

An  alternative  approach  is  provided  by  Tabarrok  [34,35]  who  formulates 
the  vibration  problem  according  to  Toupin’s  Principle  [3G]  with  the  impulse 
tensor  and  velocities  as  field  variables.  This  formulation  is  consistent  with 
the  hybrid  stress  model;  accurate  frequency  predictions  are  obtained  by 
Sakagucni  for  plate  vibration  problems  [37].  Unfortunately,  this  approach  pro¬ 
duces  neither  a  mass  matrix  nor  a  stiffness  matrix,  but  rather  a  frequency 
matrix;  the  zeros  of  its  determinant  are  the  frequencies  of  a  system.  As  a 
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result,  this  approach  is  not  applicable  to  transient  dynamic  problems. 

Recently,  efforts  have  been  made  to  develop  a  variational  principle  to 
justify  the  approach  of  interpo] ating  boundary  displacements  of  an  element  as 
a  legitimate  way  of  obtaining  a  mass  matrix  [3B,39J  so  that  the  conventional 
formulation  for  stiffness  and  mass  matrices  can  bo  safely  adopted  in  a  hybrid 
stress  finite-element  dynamic  analysis.  In  the  following,  a  brief  description 
of  the  derivation  of  a  mass  matrix  from  a  variational  principle  for  the  hybrid 
stress  model  is  given. 

Consider  a  functional  in  the  form  of  a  modified  Rcissner  principle 
[39,4t'i  for  the  free  vibration  of  a  continuum. 


j  T£(U£-Ui><M} 

a  Vn 


(3.1) 


where  Sv  refers  to  the  boundary  of  the  region  V  ,  a. .  and  u.  are  defined 
n  n  ig  1 

within  the  interior  volume  V  of  the  element  and  on  3V  ,  and  T.  and  u".  are 

n  nil 

the  boundary  traction  and  boundary  displacements,  respectively.  One  can 
recognise  that  one  of  the  equations  obtained  by  setting  =  0  is  that 
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(3.2) 


along  the  element  boundary  3V^,  where  v_.  is  the  direction  cosine  on  3V^.  Thus, 

if  I\  in  Eq.  3.1  is  replaced  by  o^v_.  and  integration  by  parts  is  carried  out 

for  the  term  1/20. .(u.  .  +  u,  .),  the  following  functional  is  obtained; 

13  1 • 3  3»i 
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Note  that  only  three  field  variables  are  involved;  0. .  and  u.  in  V  ,  and  u. 

lj  l  n  l 

on  3V  . 
n 

Assume  that  the  stress  distribution  in  the  interior  of  the  element 
may  be  written  as 

CT=P(3  (3.4a) 
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and  that  the  displacements  along  the  boundaries  of  the  element  are  assumed  to 


be  expressed  as 


Z  =  i  l 


(3.4b) 


as  before,  and  represent  the  displacement  field  u  in  the  interior  of  the 
clement  by 

H  —  hi  h  {3.5) 

where  N  is  a  "suitable"  interpolation  function  which  relates  u  to  the  nodal 
generalized  displacements  q.  Substituting  Eqs.  3.4  and  3.5  into  Eq.  3.3  re¬ 
sults  in 

TT  _  r-  f  I  Ati  i  /Z  .  /J^tn  9.  I  O T  — .  0  L  At/~  0  1 


7T„R=  i) 


where  H  and  are  the  same  as  defined  in  Eq.  2.6,  and 

B  -  ip'i’ydv' 

m  -  j/„ 

with  P’  representing  the  derivatives  of  P  (c, .  .). 

—  13,3 


The  stationary  conditions  of  TT  with  respect  to  variations  of  0  then 

mR  — 


yield 


-H(3  +  (  §-?)  f  =  0 


(3.8) 


By  solving  for  0  from  Eq.  3.8  and  substituting  into  Eq.  3.6,  one  obtains 


(3.9) 


where 


K  =  (£]  —  D)  H  (  —  P)  =  clement  stiffness  matrix 


7TI  = 


PhTf'/dV  =  "hybrid-rational"  mass  matrix 


(3.10) 


If  the  stress  function  P  is  chosen  such,  that 


1 .  i  —  0 


(3.11) 


then  L)  is  zero  and  the  element  stiffness  matrix  k  in  Eg.  3.10  reduces  to 


IS  = 


(3.17) 
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which  is  the  stiffness  matrix  of  the  hybrid  stress  model.  From  Eq.  3.‘J,  the 
following  equation  is  obtained  by  taking  6i[  =  0  since  6a  is  arbitrary; 


L  l  IS  %  +  ™  %  )  =  o 


(3.13) 


The  element  nodal  generalized  displacements  q  can  be  connected  to  the 


global  generalized  displacements  q*  by 


&  =  In"  1“ 


Note  also  that 


i  =  In*  I 


(3.14a) 


(3.14b) 


where  T  .is  a  transformation  matrix  relating  q  to  q*. 

~qq*  -a- 

Applying  Eqs.  3.14a  and  3.14b  to  Eq.  3.13,  the  following  equations  of 

dynamic  equilibrium  of  the  non- forced  complete  assembled  discretized  structure 


is  obtained; 


M  i  +K  i  =0 


(3. 13) 


where  M  and  K  represent,  respectively,  the  mass  matrix  and  the  stiffness  matrix 
for  the  complete  assembled  discretized  structure  referred  to  the  global  system. 

It  is  seen  that  the  final  governing  equations  of  free  vibration  analysis 
are  of  the  same  form  for  both  the  displacement  and  the  hybrid  stress  finite- 
element  model.  For  fre-'  vibration  analysis,  one  may  assume  that 


,Y*  2  c* 

%  =  t 


(3.16) 


where  to  is  the  natural  frequency  of  the  system.  Thus,  the  following  familiar 

omnlno  nvnh  1  ram  ypqnl  frnn  KfT  . 
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(K  -  u)  M)  2=0 


(3.17) 


3.2  Mass  Matrices  for  an  Element 


ily  resorting  to  the  modified  Reissner  principle  expressed  by  Eq.  3.1, 
the  procedure  explained  in  Subsection  3.1  shows  how  one  may  obtain  both  the 
hybrid-stress  stiffness  matrix  and  a  "hybrid  rational"  element  mass  matrix. 
Although  this  mass  matrix  is  "rational",  it  is  not  consistent  in  the  sense 
widely  understood  in  connection  with  the  finite  element  assumed  displacement 
method.  In  the  latter,  the  assumed  displacement  field  is  used  to  obtain  the 
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element  stiffness  matrix,  and  a  velocity  field  consistent  with  the  assumed 
displacement  field  is  used  to  represent  the  kinetic  energy  of  the  system  there¬ 
by  identifying  a  consistent  mass  matrix.  In  the  hybrid  stress  model  on  the 
other  hand,  the  displacement  field  utilized  for  constructing  the  kinetic  energy 
to  identify  a  rational  mass  matrix  is  not  employed  in  any  sense  in  the  de¬ 
termination  of  the  element  stiffness  matrix  given  by  Eq.  3.12. 


In  addition  to  the  "rational  element  mass  matrix  procedure"  explained 
in  Subsection  3.1,  one  may  also  construct  a  lumped  mass  matrix  by  commonly 
used  procedures,  or  from  a  special  case  of  the  present  procedure.  For  ex¬ 
ample,  i  the  displacement  field  interpolation  function  of  Eq.  3.5  were 
assumed  to  be  such  that  the  displacements  within  the  tributary  region  of  a 
particular  node  are  identical  to  the  displacements  of  that  node,  then  the  so- 
called  lumped  mass  matrix  [41]  results.  For  the  quadrilateral  element  con¬ 
sidered  in  this  study,  the  tributary  region  may  be  defined,  for  simplicity,  as 
one-fourth  of  the  total  plane-area  of  the  element  regardless  of  the  shape  of 
the  quadrilateral  element.  Then,  for  the  single-layer  element  ELEM1,  the 
following  diagonal- lumped  mass  matrix  m  is  obtained: 


(3.18) 

T  T 

where  A  =  plane  area,  p  =  density,  (q.)  =  (u.,v.,w.O  . ,0  . )  ,  i  =  1,2,  ...  4, 

1  1  1  l  xi  yi 

and  the  diagonal  terms  in  d  are  (H, ,»'.«!!,  ,H,)  with  >•  =  (h/2  dZ 

M  1113  3  1  Jh/2 

Z*Az. 

W-h/z 


and  H. 


Next,  consider  the  development  of  a  iiybrid-rational  mass  matrix  for 
ELEMl.  Since  the  element  boundary  displacements  are  assumed  to  be  linear  for 
ELEM1,  a  convenient  and  suitable  interpolation  for  fJ  would  be  a  bilinear  ex¬ 
pansion  in  terms  of  a  pair  of  transformed  coordinates  (£,n) • 
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w-  =  L  [(Ui+Zfy^Ni] 

1=1 

v  =  t  [lVi-2  9xi)Ni) 

i=» 

W  -  f_  (  vVL  N I ) 

l=l  (3.19) 


where 

N,  = 


N2  = 

N3  = 

N4  =  (/-*>>? 


(3.20) 


The  relations  between  the  original  coordinates  (x,y)  and  the  transformed  co¬ 
ordinates  (t,n)  are  (see  Fig.  23): 


X  =  f  %i  Ni 

y  =  t  &  m 


(3.21) 


Having  Eq.  3.19  (which  is  the  equivalent  of  Eq.  3.5),  the  mass  matrix  of  ELEM1 
can  be  obtained  by  substituting  Eq.  3.19  into  Eq.  3.7  and  performing  the  neces¬ 
sary  integrations.  The  resulting  mass  matrix  will  be  designated  as  m  and 
will  be  called  a  hybr id-rational  (HR)  mass  matrix  in  the  sense  that  the  inter¬ 
polated  displacements  within  the  element  conform  with  the  element  boundary 
displacements  and  is  derived  from  an  established  variational  principle. 


To  derive  lumped  cuid  HR  uiuSs  matrices  for  the  multilayer  element  I.LEMZ, 
it  is  only  necessary  to  note  that  for  each  layer  the  above-derived  mass  matrices 
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and  m 


1-IIR 


are  applicable.  However,  the  following  transformation 
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1,  2,  ...  n 


(3.22) 


Ol 

?e  =  Iez  h  <3-23> 

is  needed  to  transform  the  nodal  displacements  of  ELEMl  to  that  of  ELEMZ. 

Then  the  complete  mass  matrix  of  ELEMZ  can  be  obtained  by  summing  the  trans- 

T 

formed  lumped  or  HR  mass  matrix  T  ra  T  for  each  layer  through  the  thickness 
in  the  same  way  as  individual  element  stiffness  matrices  are  assembled  into  a 
system  matrix.  The  corresponding  mass  matrics  for  ELEMZ  will  be  designated 
as  m  and  n  ,  respectively. 

As  for  the  reduced  multilayer  clement  ELEMR,  the  corresponding  mass 

matrices  m  and  m  can  be  obtained  by  reducing  m  or  m  in  the 
^  R-L  R—HR  ^2'HR 

same  way  as  its  stiffness  matrix  is  reduced  from  that  of  ELEMZ  (see  Sub¬ 
section  2.3.1)  . 


3.3  Plate  and  Shell  Evaluation  Examples 

To  test  the  reliability  of  the  above- developed  mass  matrices,  several 
examples  have  been  solved  and  the  results  compared  with  other  known  solutions. 

A  computational  schema  making  use  of  the  Sturm  sequence  property  (42,43,44]  has 
been  adopted  to  solve  the  eigenvalue  problem  of  Eq.  3.17.  The  computer  program 
developed  is  capable  of  finding  vibration  frequencies  and  mode  shapes  for  any 
number  of  vibration  modes.  A  brief  description  of  the  solution  scheme  is  con¬ 
tained  in  the  Appendix  at  the  end  of  this  report. 


3.3.1  Single-Layer  Thin  Plates 

Both  the  lumped  mass  matrix  m  and  the  hybrid-rational  mass  matrix 
m1  UD  of  the  single-layer  element  ELEMl  are  used  in  solving  the  following  two 
problems  of  classical  thin  plate  vibrations: 

(a)  Simply-Supported  Square  Thin  Isotropic  Plate  (Table  3,  Fig.  24) 

The  dimensions  and  material  properties  of  the  plate  can  be  found  in 
Table  8,  together  with  the  numerical  results.  The  single-layer  ele¬ 
ment  ELEMl  is  used  and  the  lowest  frequency  found  by  using  n  and 

m  is  plotted  in  Fig.  24  for  various  finite-element  modelings  with 
^  I-HK 

uniform  meshes  up  to  LixO  in  a  quarter  of  the  plate;  shown  is  the  per¬ 
centage  error  for  each  of  the  present  predictions  of  the  lowest 
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frequency  in  comparison  with  the  exact  values  reported  in  Ref.  4  5- 
For  this  problem,  the  lumped  and  the  hybrid-rational  mass  matrices 
give  almost  the  same  degree  of  accuracy,  except  that  the  errors  are 
opposite  in  sign. 

(b)  Clamped  Rectangular  Thin  Isotropic  Plate  (Table  9,  Fig.  25) 

A  clamped  rectangular  plate  of  aspect  ratio  1:2  with  the  material 
properties  the  same  as  that  of  the  previous  example  has  been 
analyzed.  Again,  both  lumped  anu  HR  mass  matrices  have  been  tested; 
the  results  for  the  lowest  predicted  frequency  are  plotted  in  Fig.  25 
and  also  are  contained  in  Table  9.  Comparisons  are  made  with  the 
exact  solution  of  Ref.  45.  It  is  shown  in  Fig.  25  that  for  a  given 
number  of  degrees  of  freedom  (DOF) ,  the  HR  mass  prediction  is 
superior  to  that  obtained  by  using  the  lumped  mass  matrix. 

Since  by  using  the  present  solution  scheme,  the  computing  effort  for  a 
system  with  a  diagonal  mass  matrix  is  not  very  much  reduced  when  compared 
with  one  which  uses  a  mass  matrix  which  has  non-diagonal  terms,  in  subse¬ 
quent  calculations  only  the  HR  mass  matrices  will  be  used  so  that  better  ac¬ 
curacy  can  be  obtained. 

In  the  examples  presented  thus  far,  only  the  lowest  frequency  has  been 
computed.  To  test  the  ability  of  the  solution  scheme  for  higher  modes,  the 
following  problem  has  been  solved: 

(c)  Simply-Supported  Rectangular  Thin  Isotropic  Plate  (Table  10) 

A  simply-supported  rectangular  plate  of  the  same  dimensions  and 
material  properties  as  that  of  example  (b)  has  beer,  analyzed  for 
symmetrical  modes  of  vibration.  A  quarter  of  the  plate  was 
analyzed  by  applying  ELEM1  (with  the  HR  mass  matrix'  to  two  differ¬ 
ent  uniform  meshes  6x12  and  Hx8.  The  first  seven  frequencies  of 
symmetric  nodes  of  vibration  computed  and  the  errors  calculated  by 
comparing  with  the  exact  solution  (45]  are  shown  in  Table  10.  As 
can  be  seen  in  Table  10,  better  predictions  are  obtained  for  the 
lower  frequencies.  This  js  expected  since  the  accuracy  depends 
on  the  ability  of  the  elements  to  model  the  proper  deformed 
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shape  of  the  structure,  and  higher  nodes  involve  higher  wave  numbers 
of  deformation;  i.e.,  more  complicated  mode  shapes.  Thus,  better 
modeling  is  provided  by  the  6x12  mesh  than  by  tnc  3x8  mesh  for  the 
mode  designated  by  5-1  in  Table  10,  for  example,  because  more  ele¬ 
ments  are  put  in  the  direction  which  has  a  wave  number  equal  to  5 
(that  is,  along  the  long  side).  The  opposite  is  true  for  mode  1-3 
in  Table  10  since  nigl'.cr  wave  numbers  now  appear  in  the  other  (short 
side)  direction. 

In  summary,  the  better  solutions  of  the  two  meshes  for  the  first  five 
symmetric  modes  produce  no  error  in  predicted  frequency  larger  than  4%.  For 
higher  nodes,  apparc  itly,  more  elements  are  needed. 

3.3.2  Two-Layer  Thin  and  Thick  Plates 

Since  two-layer  laminates  are  usually  "unbalanced",  the  extensional 
mode  and  bending  mode  of  vibration  are  coupled.  The  following  two  problems 
arc  chosen  to  demonstrate  the  ability  of  che  present  finite-element  method  to 
predict  the  lowest  frequency  of  vibration  of  unbalanced  plates. 

(a)  Clamped  Cross-Ply  Square  Thin  Plate  (Table  11,  Fig,  26) 

For  thin  multilayer  plate  problems,  SLEI1R  can  be  used.  The  present 
problem  is  described  in  Table  11,  together  witli  the  results  of  the 
lowest  predicted  frequency  for  each  of  various  modeling  meshes.  Com¬ 
parisons  are  made  with  the  CPT  solution  of  Ref.  22  which  does  not 
include  transverse  siiear  deformation  effects.  The  present  finite- 
element  solution  converges  very  rapidly  and  the  small  discrepancy 
between  the  two  solutions  may  be  attributed  to  tiic  transverse  shear 
deformation  effects. 

lb)  Simply-Supported  Inf inite  Strip  (Table  12,  Fig.  27) 

In  Ref.  46,  exact  elasticity  solutions  are  obtained  for  two-layer 
cross-ply  infinite  strips  of  various  thickness- to-span  ratios. 

For  the  finite-element  analysis,  only  one  row  of  elements  arranged 
in  the  short-span  direction  is  needed.  In  fact,  only  eight  bLbMZ 
elements  are  used  along  the  half  span  because  of  symmetry.  Since  a 
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considerable  amount  of  local  (reformation  of  the  plate  cross  section 
is  expected,  the  two-layer  thickness  is  divided  into  four  sublayers. 

Even  so,  tne  finite  clement  predictions  of  the  lowest  frequency  for 
c.ic  higher  tnickness-to-span  ratio  still  are  about  20%  higher  than 
the  exact  values.  However,  as  can  be  seen  in  Fig.  27,  great  improve¬ 
ment  is  obtained  by  the  present  (ELEMZ,  HR-mass)  finite-element  solu¬ 
tion  over  (1)  the  classical  solutions  and  (2)  solutions  which  in¬ 
clude  rotary  inertia  and  in-plane  mass  inertia  but  neglect  the  local 
deformation  caused  by  transverse  shear  deforr.iation. 

3.3.3  A  Simply-Supported  Three-Layer  Square  Plate 

The  exact  vibration  analysis  of  a  three-layer  square  plate,  with  each 
layer  being  isotropic  has  been  done  in  Ref.  *57.  The  present  finite-element 
solution  uses,  ELEMZ  the  HR  mass  matrix,  and  a  4x4  mesh  in  a  quarter  of  the 
plate.  Uo  subdivision  of  each  layer  is  needed.  The  predicted  lowest  fre¬ 
quencies  show  excellent  accuracy,  as  can  be  seen  in  Table  13  for  various  combi¬ 
nations  of  density  ratios  and  shear  modulus  ratios.  Also  included  in  Table  13 
are  predictions  by  the  classical  plate  theory  which  neglects  transverse  shear 
defoliation;  the  latter,  obviously  docs  not  yield  accurate  results  even  for 
only  moderately  thick  plates  (total  thickness  =  0.1  side  length). 

3.3.4  A  Simply-Supported  Three-Layer  Cross-Ply  Cylindrical  Shell 

Numerical  solutions  for  Lota  exact  and  approximate  vibration  analyses 
of  a  three-layer  cylindrical  shell  simply  supported  at  both  ends  have  been  ob¬ 
tained  by  other  autiiors  [47,40).  The  material  properties  are  given  in  Table  14. 
In  tne-  present  investigation,  only  the  lowest  frequency  of  the  uxi symmetric  mode 
is  explored  in  order  to  restrict  the  present  computational  effort.  As  in  the 
previous  example,  only  eight  elements  in  a  half  span  are  useu.  Both  ELEMZ 
and  ELEilK  are  used  and  the  transformations  described  in  Subsections  2.3.1  and 
2.3.2  are  utilized  to  transform  the  flat-plate  nodal  displacements  of  ELEMZ 
anu  ELL MR  into  common  shell  nodal  displacements.  The  results  for  various 
thickness- to-length  ratios,  together  with  classical  (CST,  Kef.  49),  exact  (ET, 
kef.  47),  and  upj  ’'oxinate  (RET,  Kef.  4d)  solutions  arc  listed  in  Table  15. 
Excellent  accuracy  is  observed  for  ELLMZ  compared  with  the  Kef.  47  solution. 
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Even  1SLEMR,  which  neglects  the  local  warping  of  the  shell  cross  section,  yields 
fair  results  —  comparable  to  those  of  the  so-called  refined  shell  theory. 

3.4  Summary 

Both  lumped  and  hybrid-rational  mass  matrices  for  elements  ELEMl ,  BLEMR, 
and  ELEMZ  have  been  developed.  Because  of  their  better  accuracy  as  evidenced 
by  the  present  studies  of  natural  frequencies  of  vibration  for  various  plates, 
the  hybrid-rational  mass  matrices  are  chosen  for  future  use. 

Frequency  analysis  of  thin  plates  shows  that  FE  predictions  for  lower 
modes  of  vibration  arc  better  than  those  for  higher  modes,  mainly  because  the 
finite-element  discretization  employed  models  the  lower  modes  more  closely ; 
finer  meshing  is  needed  to  improve  frequency  predictions  for  the  higher 
modes  —  a  commonly  noted  fact.  The  same  phenomenon  exists  for  Raylcigh-Ritz 
approximate  analysis. 

For  very  thick  plates  (whose  thickness- to-side  ratio  is  in  the  order  of 
one)  subdivision  of  each  layer  is  needed  to  model  the  local  deformation  closely. 

Analyses  of  a  three-layer  plate  and  of  a  three-layer  cylindrical  shell 
in  the  present  study  show  excellent  accuracy  of  the  multilayer  clement  ELEMZ. 

The  present  formulation  provides  a  dependable  HR  mass  matrix  which  can 
be  confidently  applied  to  transient  dynamic  analysis. 
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SECTION  4 

THERMAL  STRESS  ANALYSIS 

Since  aerodynamic  heating  and/or  other  sources  of  heating  may  produce 
elevated  nonuniform  temperatures  in  aerospace  structures,  the  analysis  of  the 
attendant  effects  such  as  thermally -induced  distortions  and  stresses  is  of 
interest.  Accordingly  the  present  section  pertains  to  the  analysis  of  thermal 
stresses  and  deformations  in  the  context  of  applying  the  hybrid- stress  finite- 
element  model  ana  method.  This  analysis  is  developed  and  applied  to  several 
simple  thermal/structural  problems  involving  either  isotropic  or  anisotropic 
materials  in  single- layer  or  multilayer  configurations. 

It  is  assumed  that  a  steady-state  rype  of  nonuniform  temperature  dis¬ 
tribution  is  known;  the  attendant  static  thermal  stresses  and  distortions  are 
then  analyzed.  In  other  words,  the  temperature  problem  and  the  elasticity 
problem  are  considered  to  be  uncoupled. 

4.1  Thermal  Stress  Analysis  by  the  Hybrid  Stress-Model 

For  a  given  temperature  distribution,  the  thermal  stress  analysis  of 
an  elastic  body  by  the  displacement  finite-element  model  is  quite  straight¬ 
forward.  The  thermal  loading  can  be  transformed  into  a  set  of  equivalent 
nodal  forces  by  an  initial-strain  approach  [51,52].  For  a  hybrid-stress  model, 
equivalent  nodal  forces  can  also  be  obtained  by  an  initial-strain  approach  [8] . 

Consider  a  single  element.  The  complementary  energy  functional  to  be 
minimized  in  the  presence  of  a  temperature  distribution  over  the  entire  element 
is  shown  in  Ref.  8  to  be 

J  ( ^  Sijwj  &ij  +  t-oij  &ij  ~J  Ti  U-i  dA  (4.1) 

where  e  .  .  is  the  initial  strain  obtainable  by  the  following  formula  : 
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In  Eq.  4.2,  a^,  a and  are  the  coefficients  of  linear  thermal  expansion 
in  the  respective  three  principal  directions  of  the  material.  At  is  the  tem¬ 
perature  rise  from  a  stress-free  temperature  state,  and  £  and  m  are  the  plane 
direction  cosines  between  the  material  axes  and  the  element  axes.  In  short, 

Eq.  4.2  is  nothing  but  a  strain  transformation  formula  that  relates  the  initial 
strains  in  the  material  axes  to  chose  with  respect  to  the  element  axes. 

The  stress-strain  relation  in  the  presence  of  thermal  initial  strain  is 
Lij  =  °"m  StjKf  +  ioij  (4.3) 

The  inversion  of  Eq.  4.3  leads  to  the  following  expression  for  the  stress  0|_. 
which  corresponds  to  the  total  strain  : 


where 


O’ [j  —  (T ij  +  Op lj 


Oflij  ~  Cijtl 


(4.4) 


Here  0  .  is  the  stress  corresponding  to  the  initial  strain  and  0. .  is  the 
or  j  ij 

actual  stress. Since  usually  a  linear  distribution  of  a.  .  in  the  transverse 

13 

direction  can  be  assumed  for  plate  or  shell  structures,  it  follows  that  O'.  .  is 

13 

linear  in  the  transverse  direction.  The  actual  stress  c. .,  however,  may  not  be 

ij 

linear  in  the  transverse  direction,  since  o  . .  may  not  be  linear: 

013 


oy  —  frtj  ~  & ~oij 


(4.6) 


42 


lliMIUi..ll|i,.m>l.'  ,ui.i 

I  . 


I 


In  the  following,  equivalent  nodal  forces  will  be  derived  for  the  multi¬ 
layer  element  ELEMZ.  The  equivalent  nodal  forces  for  ELEMR  and  EEEM1  can  be 
obtained  directly  from  those  of  ELEMZ  by  appropriate  reductions  of  the  nodal 
degrees  of  freedom. 


For  the  stress  field,  one  may  assume  that 

£  =  P  P  +  g  it  ~  Lrb 


(4.7) 


for  each  layer,  where  P3  is  the  part  of  a! .  that  has  unknown  parameters  S, 

^"7  ,  |  * 

P  b  is  the  part  of  o! .  that  is  necessary  to  balance  the  initial  stress  P  8  , 
~t.^t  ij  f  ~0  M5 

so  that  all  together  the  actual  stress  o  satisfies  the  equilibrium  equations 
of  elasticity. 

Note  that  one  could  make  P  2  identical  to  P  8  .  Then  the  resulting 

~c~t  ,~o~o 

stresses  0__  would  be  represented  only  by  PjB,  whose  in-plane  stresses  are 
linear  in  z.  Since,  in  general,  the  actual  in-plane  stresses  may  not  be 
linear  in  z,  this  would  be  too  restrictive.  Thus  P^8t  is  not  made  equal  to 
P  8  ,  as  will  be  seen  later. 

— O'^C 

For  the  boundary  displacements ,  one  may  assume  that  t 

a  =  t  l  '  i4-ai 

as  before.  The  functional  of  Eq.  4.1  becomes,  in  matrix  notation. 


TTc  =  J  (  J  SLjic£  ftj %  +  SLJt{  fcj  UV-j-Ti  UL  M 

=  j~(  P  £  +  ^  ~  ftR  Pg  +  f}Pi  -  J  dV 

=  j  (5-(5Tpr5P(3+  prp'sptPtt±  fiXs 
‘  j  C  I  +  L±^X  LpJA 


(4.0) 
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In  Eg.  4.9,  the  matrices  JR,  R^_,  and  R^  are  obtained  by  specifying  boundary 

coordinates  in  the  function  P,  P  ,  and  P  ,  respectively,  of  Eq.  4.7. 

—  ~t  — o 

By  expressing 


{4.10b) 


H  =  J  pJS  P  dV  {4.10a) 

<3,  =  j  RTL  dA 

as  before,  and  defining 

Ht  =  J  PTS  Pf 

Olt^i^LdA  {4*10b) 

:$0  =  Jg°LdA 

the  functional  in  Eq.  4.9  becomes 

\  =5  PTHp  -P’SJ  fttt  +  gG.f 

C*  -V  A/  A/  — v  -Or  /V  -v  -V-  ^  'V  (411) 

In  Eq.  4.11,  the  constant  terms  /$.TP.TSP  3  dV  and  /B  TP  rsP  3  dv  are 
-  1  ~  t  -»t  ~^t*t  ~o  ••  o  ~~o*o 

dropped,  since  they  have  no  contribution  upon  taking  ,r.riations  of  7i  . 

c 

Since  a  multilayer  element  is  being  considered,  the  interlayer  equilibrium 
of  stresses  will  lead  to  a  set  of  constraint  on  6  (see  Subsection  2.3.1): 


(4.11) 


A  (3=  3 


(4.12) 


Upon  introducing  the  Lagrange  multiplier  A,  a  new  functional  is  formed  by 
combining  Eq.  4.11  and  Eq.  4.12: 


7T  *  -  Jlc  +  ( A  (3  -  8  ) 


(4.13) 


Taking  the  variation  of  n*  with  respect  to  3,  one  obtains 

c  ^ 


H  |3  -  0]  ^  +  Ht  (3.  -f  At>  =  q 


(5  =  H',(6$-WtPt -ATA) 

.  -v  — w  X  ■— W  -V  A/  -V 


(4.14) 


Substituting  6  of  Eq.  4.14  into  Eq.  4.12,  one  obtains 
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AH"'(<j  !'Ht£+-ATA)  =  B 

or 

A  =  (  AHV/'(AH;'5_|-AH‘,Ht£t-6)  (4.15) 

Substituting  3  and  A  back  into  Eq.  4.13,  the  new  functional  it’  now  contains 
1  c 

only  q: 


He  =  -{ 

-  $VwV (A HVj'  +  p,r6HH,  fit ) 

+  fCiJ  ft, 

/V  /V-  -w 

(4.16) 


Again,  a  constant  term  in  TP  has  been  dropped.  Finally,  by  taking  the  varia¬ 
tion  of  iP  with  respect  to  q,  one  obtains 


=  -e,'TH*,ATUH-A")_'(B+AH-,H4ft)-(<5t',fe-  SVHt( \)+<£(S0 

^  v  -V  -V  -v  -V  -v 

(4.17) 

Thus,  the  element  stiffness  matrix  Jk  of  Eq.  2.28  remains  valid  and  an  equiva¬ 
lent  element  nodal  thermal  force  vector: 


f  =  CiT H-'«t  ft ~cj  HAr(A  H*At)'(&+AH\ ft )  +  0,1  ft  -  Qt  ft 

'V  Ay  /v  aj  ^  ^  ~  ~  -v  -v  -u  •v'  'v  ^ 


(4.18) 


is  obtained. 

4.2  Modeling  of  a  Temperature  Distribution 

The  equivalent  tnermal  loading  vector  derived  in  Subsection  4.1  depends 

on  the  temperature  distribution  of  an  element.  In  theory,  any  given  function 

of  temperature  distribution  can  be  used  in  Eq.  4.2  to  yield  a  vector  of  initial 

strain  and  eventually  the  vector  P  l3  in  Eq.  4.7.  That,  however,  would  intro- 

o  o 

duce  difficulties  in  evaluating  Eq.  4.10b.  Since  the  basic  stress  distribution 
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represented  by  £8  in  Eq.  4.7  is  only  linear  in  x  and  y  (see  Subsection  2.3.2), 
higher-order  descriptions  of  temperature  for  the  equivalent  nodal  force  calcu¬ 
lation  may  be  essentially  futile  unless  higher-order  stress-distribution 
functions  are  used  in  developing  the  hybrid  stress  finite-element  properties. 
Therefore,  the  following  simplification  in  temperature-distribution  representa¬ 
tion  is  adopted  in  this  study.  A  given  arbitrary  temperature  distribution 
AT^(x,y,z)  for  any  one  layer  of  an  element  is  modeled  by  a  simple  distribution 
AT{x,y,z) : 


AT(*.y,zj  =  T,  +  Tz*  +T3y  +  *(  T4  +  T*  *+•  T&yJ 
+  Z2  (  Tj  ■+■  Te:x  +  T9  y) 


(4.19) 


which  is  parabolic  in  z  and  linear  in  x  and  y.  The  coefficients  T  ,  T^,  ...  Tg 
are  determined  by  minimizing  the  square  errors  of  surface-fitting  AT  at  the 
12  corner  points  of  that  layer  of  an  element  (Fig.  28a)  . 

Equation  4.19  can  also  be  expressed  in  terms  of  temperatures  in  the 
bottom  (T^) ,  middle  (T^) ,  and  top  (Tt>  surface  of  that  layer  (see  Fig.  28b); 

T  +  Tz*  +T3y  =  Tm,  +  Tmz?  +■  rm3  y 

T4  +  T5*  +  Tfcy  =  1  CrTt,-Tbt)+rTt2-T^)*+fTt3-TM)yJ 

T7  +  T8y  +T9  y=  ^  r  ( Ttl  -Thl  -  2 Tn, )  +  ( ru -r^-2  r** )  x 

+  (  T Ty  ~2Tmi )y) 


(4.20) 


where  h  is  the  thickness  of  that  particular  layer. 


Thus,  the  minimization  of  errors  in  AT(x,y,z)  can  be  carried  out  for  the 
three  surfaces  separately.  Denote  by  T  the  actual  temperatures  at  the  four 
corner  points  for  any  one  of  the  three  surfaces  calculated  from  the  original 
AT^(x,y,z) ,  and  denote  by  T\  the  approximate  values  calculated  from  Eq.  4.19 
and  Eq.  4.20.  Then  the  square  error  E^_  of  the  temperature  description  for  that 
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surface  is 


Er=  &(VT.) 


1  =  1 


and 


Tt  Ts,  +  T52  *L  +  t$3  y(‘ 


s  =  b,  or  m,  or  t 


(4.21) 


(4.22) 


where  x.  and  y.  are  the  corner  coordinates  of  the  element. 
1  1 

To  minimize  E  : 
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Y?  =  t  -IT)  =  O 

0  Is}  i-J 


(4.23) 


After  rearranging  the  terms  in  Eq.  4.23,  and  using  Eq.  4.22,  the  following 


equations  are  obtained: 
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(4.24) 


Then  T  , ,  T  ,  T  ,  can  be  solved  easily  in  terms  of  temperatures  at  corner 
si  s2  s3 

points  T0^>  and  the  modeling  of  the  temperature  distribution  AT^(x,y,z)  by 
AT(x,y,z)  is  completed  upon  substituting  the  expressions  in  Eq.  4.20  into 
Eq.  4.19. 
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4.3  Stress  Assumptions  in  an  Element 


The  actual  stress  o  consists  of  three  parts:  P0,  P  b. ,  and  P  0  as 

—  ~-t~t  ~-o~o 

can  be  seen  in  Eq.  4  7,  For  each  layer  of  a  multilayer  element,  P0  is  the 

same  as  that  of  ELEM1,  and  P  0  can  be  calculated  from 
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(4.25) 

where  the  matrix  C  is  the  elastic  matrix  and  can  be  obtained  by  inverting  the 

compliance  matrix  S  of  Eq.  2.14;  the  initial  strain  vector  is  obtained  from 

Eq.  4.2.  Since  P0  satisfies  the  three  equilibrium  equations  of  elasticity 

identically,  it  is  now  necessary  to  derive  P0_  such  that  ( P  0  -  P  0  )  will 

~t^t  t~t  ~o~o 

satisfy  the  equilibrium  equations  identically.  Substituting  P  0  of  Eq.  4.25 

o-o 

into  the  following  equilibrium  equations 

(  +  4  (-G-tyzl.Z  —  o 

(4.26) 

and  integrating,  one  finds  the  following  for  £^0^: 
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(4.2?) 

Thus,  the  stresses  of  any  layer  of  an  element  consist  of  £J3  of  ELEM1, 

P  3  of  Eq.  4.23,  and  P  3.  of  Eq.  4.27. 

-O'-o  *  ''-t.M: 

The  interlayer  equilibrium  conditions  impose  constraints  on  0^,  c  * 
and  G  (see  Eq.  2.30).  The  equilibrium  of  G^  is  automatically  satisfied, 
since  G  of  P3  is  zero  and  G  =  0  +  0  -G  =0byEq.4.27.  The  equi- 

2  2  tZ  OZ 

librium  of  0  and  G  leads  to 
xz  yz 

C  ( <T*2 , 6y2 )  of  P(3  +((^2.^)0^ at  the  top  of  the  ith  layer 

(4.28) 

=  (  (trx2,(Ty2)of  Pfl  +  c<rx2 .  <j>2  at  the  bottom  of  the  (i+l)^ 

layer 

since  P  3  does  not  contain  nonzero  transverse  shear  stresses.  In  Eq.  4.28, 
-o-'o 

the  part  of  Pg  yields  the  (same)  coefficient  matrix  £  of  3  in  Eq.  21  of 

ELEMZ,  and  the  part  of  P  3  yields  constant  terms  which  constitute  the  vector 

t^t 

3  in  the  constraint  equation 

A  {3  =  6  (4.29) 

Thus,  all  of  the  ingredients  needed  to  calculate  the  equivalent  nodal  forces 
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are  now  available.  In  summary,  the  equivalent  nodal  forces  can  be  obtained 

from  Eq.  4.18,  in  which  the  matrices  G^^and  are  the  same  as  defined  for 

ELEMZ  before,  the  vectors  H  6  and  G  3  are  obtained  by  using  Eq.  4.10  and 

T  ^ 

Eq.  4.27,  the  vector  G  3  is  obtained  from  Eq.  4.10  and  Eq.  4.25,  and  the 

~-o  ~o 

vector  £  is  obtained  from  Eqs.  4.27  and  4.28. 


4.4  Plate  and  Shell  Evaluation  Examples 

4.4.1  In-Plane  Expansion  of  a  Free  Single-Layer  Thin 
Isotropic  Rectangular  Plate 

To  test  the  accuracy  of  the  equivalent  thermal  loading  representation, 
a  problem  for  which  both  experimental  and  approximate-theoretical  [53]  thermal 
stress  data  are  available  is  chosen  as  a  test  problem.  It  is  a  free  single¬ 
layer  thin  rectangular  isotropic  plate  heated  linearly  in  the  plate  plane  and 
unif'orml  •  across  the  plate  thickness  direction  as  illustrated  in  Fig.  29;  other 
pertinent  data  for  this  problem  are  given  in  Table  16.  Because  of  the  sym¬ 
metry  of  the  plate  and  the  temperature  distribution,  only  a  quarter  of  the 
plate  is  used  in  the  analysis  and  only  the  in-p^ane  displacements  and  stresses 
need  to  be  considered. 


The  in-plane  displacement  pattern  calculated  by  the  present  finite- 
element  method  using  a  5x4  mesh  is  plotted  in  Fig.  30.  In  Figs.  31a,  31b,  and 
31c,  the  in-plane  stresses  o^,  0^,  and  0  ,  respectively,  are  plotted  against 

both  experimental  and  theoretical  results  of  Ref.  53.  It  is  seen  that  the 
direct  stresses  (0^  and  o^)  calculated  by  using  the  present  finite-element 
method  agree  very  well  with  both  results;  the  shear  stress  (o^J  prediction 
is  also  very  good,  except  near  the  edges,  because  the  free  traction  edge  con¬ 
dition  is  satisfied  only  in  ^n  average  sense  in  the  present  finite  element 


modeling.  It  should  be  mentioned  that  the  theoretical  results  of  Ref.  53  are 
only  approximate;  the  direct  stress  0^  is  first  assumed  according  to  a  dis¬ 
tribution  which  is  exact  for  an  infinite  strip  under  thermal  loading  and  the 
other  stresses  are  then  determined  according  to  the  Principle  of  Complementary 
Energy.  Also,  as  mentioned  in  Ref.  53,  the  temperature  control  to  give  the 
temperature  distribution  of  Fig.  29  for  the  experiment  is  not  perfect.  There¬ 
fore,  some  slight  discrepancies  between  the  theoretical  and  experimental 
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results  are  expected.  However,  the  present  prediction  agrees  very  well  with 
both  of  the  comparison  results.  This  demonstrates  the  accuracy  of  the  present 
equivalent  thermal  loading  calculation  and  the  associated  FE  analysis. 

4.4.2  Bending  of  a  Clamped  Thin  Isotropic  Rectangular  Plate 
Subjected  to  Monuniform  Temperature 

The  dimensions  and  properties  of  the  plate  under  consideration  are 
listed  in  Table  17.  For  a  thin  rectangular  isotropic  plate  under  thermal 
loading  and  having  all  four  sides  ideally  clamped,  the  thermal  stress  problems 
can  be  solved  by  recognizing  that  the  governing  equations  on  the  transverse 
displacements  are  identical  to  those  of  the  same  plate  under  distributed 
static  lateral  loading.  In  Ref.  54,  it  is  shown  that  this  equivalent  static 
load,  denoted  by  p*  is 

P*=-  7“  (VJMt)  (4*30) 

where 

Mt=  a^l/z7ZjZ  (4.31) 

where  v  is  the  Poisson  ratio,  E  is  Young's  modulus,  a  is  the  linear  thermal 
expansion  coefficient,  and  T  is  the  temperature  distribution.  In  the  present 
example,  a  temperature  distribution  constant  in  x,  parabolic  in  y,  and  linear 
in  z  is  assumed  (see  Fig.  32a).  Thus,  the  results  calculated  by  the  present 
hybrid-stress  FE  approach  can  be  compared  with  that  of  the  classical  solution 
(18]  for  a  clamped  rectangular  plate  under  uniform  loading.  The  displacements 
along  the  central  lines  of  the  plate,  obtained  by  the  present  method  with  an 
8x3  uniform  mesh  in  a  quarter  of  the  plate  are  plotted  in  Figs.  32b  and  32c, 
together  with  the  central  deflection  predicted  by  the  classical  solution  for 
static  bending  (18] .  The  present  prediction  for  the  central  deflection  is 
almost  identical  to  the  classical  exact  solution  of  Ref.  18. 

4.4.3  Simply-Supported,  Infinitely-Long,  Two-Layer,  Cross-Ply 
(0o/90o)  Thin  Strip  Subjected  to  a  Uniform  Temperature 
Distribution 

The  material  properties  and  dimensions  of  the  thin  flat  strip  are  de¬ 
scribed  in  Table  18;  the  pertinent  coordinates  and  geometry  are  depicted  in 
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Fig.  33a.  Since  this  laminated  plate  is  not  balanced  in  the  sense  that  the 
expansion  properties  are  not  symmetric  in  z,  even  a  uniform  temperature  will 
cause  a  transverse  displacement  as  can  be  seen  in  Fig.  33b.  Also,  thermal 
stresses  are  induced.  Plotted  in  Figs.  33c  through  33f  are  the  in-plane 
stresses  and  interlaminar  shear  stresses.  No  other  solutions  are  available 
at  present  for  this  problem;  therefore,  no  comparisons  can  be  made.  From 
Figs.  33c  and  33d  it  is  seen  that  the  maximum  occurs  near  the  simply- 

supported  edges  (x  =  +_  8  in. )  in  the  lower  layer  because  the  lower  layer  has 
a  larger  elastic  constant  in  che  x-direction.  The  maximum  interlaminar 
shear  stress  (Fig.  33e)  occurs  at  the  supported  edge  and  the  stress 

is  almost  constant  along  the  x-direction. 

4.4.4  Simply-Supported,  Three-Layer,  Cross-Ply  (0°/90l>/0<> )  , 

Thick,  Cylindrical  Shell  Subjected  to  a  Uniform 
Temperature  Distribution 

A  uniform  temperature  is  applied  to  a  three-layer  cylindrical  shell  of 
finite  length,  simply  supported  at  each  end.  The  material  properties  and  the 
dimensions  of  the  shell  are  described  in  Table  19.  Since  this  is  an  axisym- 
metric  problem,  only  one  strip  of  the  shell  is  needed  in  the  finite-element 
modeling.  In  fact,  only  half  a  strip  is  used  because  of  symmetry  in  the  axial 
direction  also.  The  calculated  radial  displacements  resulting  from  the  applied 
temperature  are  plotted  in  Fig.  34a.  It  is  interesting  to  note  that  the  maxi¬ 
mum  radial  deflection  does  not  occur  at  the  cei  tral  section  but  at  a  place  near 
the  center  section.  However,  this  phenomenon  is  a  result  of  the  particular 
combination  of  the  length,  thickness,  and  radius  ratios  used  in  this  example 
and  should  not  be  regarded  as  a  general  result.  In  Fig.  34b,  the  axial  dis¬ 
placements  of  six  representative  sections  are  plotted.  The  distortion  from 
the  original  plane  section  is  most  severe  at  the  supported  edge  but  dies  out 
gradually  toward  the  central  section.  A  reverse  of  the  direction  of  rotation 
of  tne  cross  section  occurs  near  the  central  section  (x=7) .  This  is  consistent 
with  the  small  decrease  of  radial  displacement  near  the  central  section 
(Fig.  34a). 


52 


4.5  Summary 

The  thermal-stress  analysis  by  the  present  hybrid-stress  finite- 
element  method  is  shown  to  be  equivalent  to  a  static  analysis  Cas  is  well 
known  to  be  true  also  for  other  types  of  thermal  analyses) .  The  equivalent 
thermal  loading  vector  for  a  hybrid  stress  model  is  derived  using  an  initial- 
strain  approach.  Arbitrary  temperature  distributions  on  any  element  are  ap¬ 
proximated  by  simple  functions:  linear  in  x  and  y,  and  parabolic  in  z.  The 
accuracy  of  the  present  method  is  demonstrated  in  two  simple  problems:  (1) 
in-plane  plate  expansion  and  (2)  thermal-induced  bending  of  a  plate;  good 
agreement  is  observed  for  both  displacement  and  stress  predictions.  Two 
other  problems,  a  two-layer  plate  and  a  three-layer  cylindrical  shell,  each 
under  a  uniform  temperature  distribution,  are  also  solved  and  the  results  are 
discussed. 


SECTION 


SUMMARY  AND  CONCLUSIONS 


5.1  Summary 

This  study  has  been  devoted  to  the  development  and  evaluation  of  hybrid 
stress  flat-plate  finite  elements  for  use  in  analyzing  laminated  plate  and 
shell  structures.  Transverse  shear  deformation  effects  as  well  as  in-plane 
and  rotary  inertia  for  vibration  analysis  are  included.  The  plate  and  shell 
structures  under  consideration  range  from  very  thin  to  very  thick.  According 
to  the  loading  conditions  applied  to  the  structure,  this  study  can  be  divided 
into  three  categories:  static  analysis,  vibration  analysis,  and  thermal  stress 
analysis. 

These  three  types  of  loading  conditions  are  discussed  since  they  repre¬ 
sent  trie  basic  loading  environment  that  many  structures,  particularly  aero¬ 
nautical  and  aerospace  structures,  may  encounter.  Tnese  basic  analyses  serve 
as  the  foundation  for  future  studies  on  transient  dynamic  and  other  compli¬ 
cated  loading  conditions.  Also,  they  are  necessary  aspects  in  the  evaluation 
of  the  properties  of  certain  finite  elements  developed  so  that  the  reliability 
and  limitations  of  these  finite  elements  can  be  assessed  before  these  finite 
elements  are  applied  to  practical  structural  analysis  and  design. 

In  static  finite-element  analysis,  the  basic  task  is  to  develop  reliable 
stiffness  matrices  that  converge  rapidly  to  correct  solutions,  with  a  "minimum" 
number  cf  degrees  of  freedom.  To  achieve  this  goal,  studies  were  made  first  on 
single- layer  elements.  An  outline  of  the  derivation  of  single-layer  flat- 
plate  elements  by  the  hybrid  stress  mouei  is  given  and  a  LoLai  of  nine  c’ements 
witn  different  stress  and  boundary  displacement  assumptions  were  tested  on  a 
very  thin  flat  plate  as  well  as  on  a  very  thick  flat  plate.  Among  these  nine 
finite  elements,  a  quadrilateral  element  designated  as  KLEM1  {with  linear  boundary 
displacements  and  linear  ir.-plane  stress  assumptions)  is  enosen  for  use  because 
of  its  fast  convergence,  excellent  accuracy,  and  simplicity. 

Then  trie  derivation  or  the  stiffness  matrix  for  a  multilayer  element 


was  carried  oat  by  extending  the  basic  sinyle-layer  flat-plate  element  KLEKl 
to  a  multilayer  element  (designated  by  ELEMZ)  which  has  independent  boundary 
displacements  that  can  model  the  local  deformatio’  of  each  layer  of  a  plate 
or  shell  closely.  For  relatively  thin  multilayer  plates  or  shells,  an  ele¬ 
ment  designated  by  ELEMR  is  obtained  by  reducing  the  nodal  degrees-of- freedom 
of  hum. 

Two  transformation  formulas  are  given  for  sinyle-layer  and  multilayer 
fiat-plate  elements  so  that  their  appl-^ation  can  be  broadened  to  include 
curved-snell  problems. 

Various  static-plate  and  curved-shell  problems,  including  both  thin  and 
thick,  single-layer ,  two-layer,  and  three- layer  problems,  have  been  solved 
using  the  three  elements  developed  in  this  study,  and  comparisons  have  been 
made  with  other  existing  solutions.  It  is  founa  that  tnese  elements  are  indeed 
very  effective  in  predicting  both  displacements  and  stresses  of  plate  and  shell 
structures  under  static  loading. 

For  vibration  analyses  using  the  hybrid-stress  model,  a  brief  discussion 
of  the  variational  basis  of  the  derivation  of  a  mass  matrix  is  given,  followed 
.by  the  development  of  lumped  and  hybrid -rational  mass  matrices  for  the  three 
elements  mentioned  earlier.  Based  on  the  results  of  test  problems,  the  less 
accurate  lumped-mass  approach  was  discarded,  and  the  superior  hybrid-rational 
mass  approach  was  adopted  for  subsequent  use.  Several  example  problems  were 
solved.  Frequency  predictions  were  compared  with  other  existing  solutions, 
and  good  accuracy  of  the  present  predictions  is  observed  in  general. 

The  steady-state  thermal  stress  analysis  carried  out  by  the  hybrid-stress 
finite-element  model  is  characterized  by  'i.e  development  of  an  equivalent  therm., 1 
loading  vector;  tne  analysis  is  otherwise  1 imilar  to  a  static  analysis.  The 
temperature  distribution  wi-„.iin  an  element  is  modeled  for  convenience  in  the 
present  study  by  a  linear  approximation  in  the  plane  surface  and  by  a  parabolic 
fit  in  the  transverse  direction.  Then  initial  strains  and  their  correspond’  ig 
stresses  are  computed  based  on  these  temperature  distributions;  the  total  stress 
assumptions  are  such  that  they  satisfy  the  equi librium  conditions  of  elasticity. 
Tne  equivalent  thermal  loading  vector  is  cotaii.ed  through  the  application  of  a 
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variational  principle  that  includes  the  temperature  effects  as  initial  strains. 

Problems  involving  in-plane  expansion  and  transverse  bending  induced  by 
temperature  were  solved.  Because  of  the  limited  number  of  existing  numerical 
solutions  on  thermal  stress  problems,  only  two  of  the  example  problems  were 
compared  to  other  solutions;  they  show  very  good  accuracy  for  displacement  and 
stress  predictions. 

In  summary,  three  quadrilateral  flat-plate  elements:  ELEM1,  ELEMR,  and 
ELEI1Z  were  developed  and  verified  to  provide  reasonably  accurate  and  reliable 
predictions  for  static,  vibration,  and  thermal  stress  problems.  They  all  have 
similar  assumptions  on  stresses  and  boundary  displacements  but  differ  in  their 
applicability  to  (1)  thin  or  thick  and  (2)  single-layer  or  multilayer  plate 
and/or  shell  problems. 

5.2  Conclusions 

Based  upon  the  results  of  various  numerical  verifications  included  in 
thi^  study,  the  following  conclusions  may  be  stated: 

(a)  The  present  finite-element  method  is  reasonably  efficient 
and  accurate  for  predicting  stress  and  displacement  re¬ 
sponses  to  static  or  thermal  loadings  and  natural  fre¬ 
quencies  for  free  vibrations  of  plate  and  shell  structures. 

(b)  For  single-layer  problems,  ELEMl  is  recommended  for  use. 

However,  if  the  thickness  becomes  large  (for  thickness- 
to-side  ratio  h/l  greater  than  about  0-1),  it  is  advisable 
to  subdivide  the  single-layer  into  sublayers  and  use  ELEMZ 
to  model  the  local  deformation  more  closely. 

(c)  For  relatively  thin  multilayer  plates  or  shells  with  h/l 
or  h/R  (R  =  radius)  smaller  than  0.1,  the  ELEflR  element 

can  be  used.  It  lias  been  pointed  out  in  Subsection  2.7  that 
the  element  reported  in  Ref.  29  is  a  better  choice  for  this 
range  of  application. 
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1  : 

(d)  For  (1)  relatively  thick  plates  or  shells  or  (2)  laminates 
that  have  _elativciy  low  shear  moduli  and/or  a  !high  ratio  of 
elastic  constants  in  different  directions,  the  ELEMZ  multilayer 
element  should  be  used  so  that  local  deformations  and  stresses 
arc  predicted  accurately. 

(e)  Predictions  or  the  lowest  natural  frequency  are  usually 
accurate  even  when  a  small  number  of  elements  is  used', 

v>'  .  i 

but  more  elements  are  needed  for  suitably  accurate  predic¬ 
tions  of  the  higher  frequencies.  * 

(f)  Comparisons  of  predictions  from  using  two  assumed-displace¬ 
ment  finite  elements  for  thick  laminated  plates  on  static 

problems  snowed  that  the  present  hybrid-stress  element 

.  '  '  '  1 

ELEMZ  provides  more  accurate  preaictidns  for  a  given  number 

of  unknowns  (and  hence  is  less  expansive) . 1 

(g)  Thermal  effects  are  accurately  accpunted  for  through  the 
computation  of  an  equivalent  thermal  loading  vector  in 
the  present  hybrid-stress  formulation.  The  limitations 

on  tiie  present  finite-element  methods  lie  in  the  fact  that 

the  normal  stress  a  is  assumed  to  be  zero  and  the  transverse 
z 

displacement  w  is  assumed  to  be  constant  throughout  the 

whole  thickness  of  an  element.  While  errors  in  the  o  or  the 

z 

w  prediction  are  usually  negligible,  these  restrictive  as¬ 
sumptions  do  have  some  effect  on  the  other  stress  predictions. 

The  present  approach  always  yields  symmetric  results  witli  re¬ 
spect  to  the  middle  surface  of  an  element  if  the  laminate  is 
symmetric;  however,  the  actual  results  may  involve  some  devi¬ 
ation  from  symmetry  because  the  load  may  act  on  only  one  sur¬ 
face  of  the  laminate.  Fortunately,  such  deviation  is  usually 
negligible  except  for  very  thick  plates,  for  example,  for  h/l 
greater  than  about  0.25  [ 2 J . 

In  general,  tUu  stiffness  matrices,  m.,ss  matrices,  and  the  equivalent 
thermal  loading  vectors  aeve loped  in  this  study  are  very  effective.  They  can 
be  used  for  detailed  displacement  and  stress  analysis.  Also  they  constitute  a 
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sound  foundation  for  transient  dynamic  analysis. 

5.3  Suggestions  for  Future  Research 

A  direct  and  useful  extension  of  the  present  work  would  be  the  study  of 
the  response  of  layered  structures  under  transient  dynamic  loadings.  Since 
the  basic  stiffness  matrix  and  the  mass  matrix  are  already  developed,  in  theory 
it  is  necessary  only  to  select  an  efficient  timewise  numerical  integration 
scheme  [55]  to  solve  the  transient  equations  of  motion  if  damping  can  be  neg¬ 
lected.  Otherwise,  a  damping  matrix  will  have  to  be  developed. 

The  geometric  stiffness  matrix  needed  in  a  buckling  analysis  could  be 
developed  in  a  manner  similar  to  the  treatment  used  in  the  derivation  of  a 
mass  matrix  by  the  hybrid  stress  model.  In  fact,  the  derivation  of  a  geometric 
stiffness  matrix  based  on  a  mixed  variational  principle  has  already  been  reported 
by  Allman  [56];  even  earlier  Lundgren  [57]  obtained  a  geometric  stiffness  matrix 
for  laminated  plates  by  independently- assumed  displacement  functions.  Thus,  the 
extension  of  the  present  work  to  a  linear  buckling  analysis  including  thermal 
buckling  seems  to  be  quite  feasible. 

Another  area  of  interest  is  the  nonlinear  analysis  of  plates  and  shells. 
Both  geometric  and  material  nonlinearities  could  be  considered.  For  geo¬ 
metrically  nonlinear  problems,  the  incremental  approach  [5B]  seems  to  be  the 
most  direct  solution.  However,  some  modifications  of  the  present  finite -element 
models  would  be  needed  since  the  present  quadrilateral  flat-plate  element  will 
not  fit  the  deformed  shape  of  a  plate  or  a  shell.  For  elastic-plastic-type 
material  properties,  some  progress  has  been  made  recently  [59]  for  hybrid 
stress  finite-element  models.  The  inclusion  of  elastic-plastic  behavior  in 
this  hybrid-stress  finite-element  context  may  turn  out  not  to  be  as  difficult 
as  the  treatment  of  geometrically  nonlinear  problems.  These  nonlinear  effects, 
of  course,  have  been  accommodated  in  assumed-displacement  finite-element,  finite- 
difference,  and  other  methods  for  isotropic  single-layer  and/or  multilayer 
structures. 
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Q,Q  =  Respectively,  complete  and  incomplete  quadratic  functions  of  x  and  y 
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TABLE  2 


DATA  FOR  SINGLE-LAYER  ISOTROPIC  PLATES* 


(a)  Simply-Supported  Thin  Plate 

Loading:  Concentrate  Central  Loadiny 

Aspect  Ratio  =  1:1 

Thickness- to-siue  Ratio  =  1:100 

Boundary  Condition:  Simply-Supported,  no  tangential  rotation 
Finite  Element  Breakdown:  Uniform  mesh  in  a  quarter 

(b)  Clamped  Square  Thin  Plate 

Loading:  Concentrate  Central  Loading 
Aspect  Ratio  ="  1:1 
Thickness-to-side  Ratio  =  1:100 
Boundary  Condition:  Clamped 

Finite  Element  Breakdown:  Uniform  mesh  in  a  quarter 


(c)  Clamped  Rectangular  Plate 

Loading:  Uniform  Load 
Aspect  Ratio  =  1:2 

Thickness-to-snort-side  Ratio  =  1:100 
Boundary  Condition:  Clamped 

Finite  Element  Breakdown:  Uniform  mesh  in  a  quarter  with 

element  aspect  ratio  1:1  and  1:2 

(d)  Skewed  Simply-Supported  Thin  Plate 

Loadiny:  Uniform  Load,  p(psi) 

Aspect  Ratio  =  1:1 
Thickness-to-side  Ratio  =  1:100 
Length  of  each  side  =  a 
Boundary  Condition  =  Simply-Supported, 

Finite  Element  Breakdown:  Uniform  mesh 

(e)  Simply-Supported  Square  Thick  Plate 

Loading:  Uniform 

Aspect  Ratio  =  1:1 

Thickness-to-side  Ratio  =  1:4 

Boundary  Condition;  Simply-Supported,  no  tangential  rotation 

Finite-Element  Breakdown:  Uniform  mesh  in  a  quarter 

*7  3/2 

L  =  10  usi,  v  =  0.3  for  all  cases;  0  =  En/[l2(l  -  V  )],  li  =  plate  thickness 
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tangential  rotation  allowed 
for  the  whole  plate 


TABLE  3 


DATA  FOR  TWO- LAYER  PLATES 


(a)  Clamped  Cross-Ply  Rectangular  Thin  Plate 
Loading:  Uniform 
Aspect  Ratio  =  1:2 
Tliickness-to-side  Ratio  =  1:100 
Boundary  Condition:  Clamped 
Finite  Element  Breakdown:  Uniform  mesh  in  a  quarter,  square  elements 
Material:  =  40,  G12/E2  =  0.5,  v  =  0.25 

Orientation:  0°/90° 


(b)  Clamped  Cross-Ply  Square  Thin  Plate 

Same  as  in  (a),  except  Aspect  Ratio  =  1:1 

(c)  Clamped  Angle-Ply  Square  Thin  Plate 

Same  as  in  (b) ,  except  fiber  orientation  =  +  45° 

(d)  Simply-Supported  Angle-Ply  Square  Thick  Plate 

Loading:  q^  sin(irx/a)  sin(Tiy/a) 

Aspect  Ratio  =  1:1 
Thickness-to-side  Ratio  =  1:6 

Boundary  Condition:  Simply-Supported,  free  in  tangential  direction, 

no  normal  displacement 

Finite  Element  Breakdown:  Uniform  mesh  for  the  whole  plate 

Material:  E,/E„  =  40,  G  /E  =  O.b,  G^VE^  =  0.5,  V,  =  0.25 
1  2  12  2  23  2  12 

Orientation:  6  varies 
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TABLE  4 


Orientation:  Sane  as  in  (a) 

(1)  Aspect  Ratio  =  3:  b  =  3a 
Tnickness-to-short-side  Ratio  h/a  =  1/10 

(2)  Aspect  Ratio:  a/b  =  1 


Thickness-to-side  Ratio  h/a  =  1/4 

Nondimensionalized  Quantities 

u  *  u(a/2,  0,z)  ET/{q  a3/h2) 

_ox  «  a  (0,0,z)/{q  a2/h^) 

0xy  =  0*y(a/2,  a/2?  z)/(qoa  /ti  ) 

*xz  =  °xzU/2'  °'  z)/tV/h) 


od 


TABLE  6 


DATA  FOR  THIN  CYLINDRICAL  SHELLS 


(a.  Pinched  Cylindrical  Shell 

Loading:  Pinch  Load,  P  (lbs) 

Radius:  a  =  10  in. 

Length  l  =  31.42  in. 

Thickness:  h  *  0.1  in. 

Boundary  Conditions:  Simply-Supported,  no  tangential  rotation 
Finite  Element  Breakdown:  Uniform  mesh  in  one-eighth 
Material:  Isotropic,  E  =  1.0  psi,  V  =  0.28 


(b)  Ring-Loaded  Cylindrical  Shell 

Loading:  Axisymmetric  Ring  Load  of  Q  =  1.0  lb/in 
Radius:  a  =*  10  in. 

Length:  i  ^  31.42  in. 

Thickness:  h  =  0.1  in. 

Boundary  Conditions:  Free  edge 
Finite  Element  Breakdown:  Uniform  mesh  in  one  strip 
Material:  Isotropic,  E  =  1.0  psi,  V  =  0.28 
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TABLE  7 


DATA  FOR  CONICAL  SHELLS 

(a)  Isotropic  Single-Layer  Conical  Shell 

Loading:  Axisymmetric  King  Loads  near  free  end f  Q 
Radius  at  smaller  end  =  10  in.  1 

Radius  at  larger  end  =  15  in. 

Length  of  meridian  i  10  in 
Semi-angle  of  cone  =  30° 

Thickness  =  0.025  in.,  0.1  in.,  0.5  in.,  1.0  in.  .  . 

Boundary  Condition:  Clamped  at  smaller  end,  free  at  larger  end 

7 

Material  =  Isotropic,  E  =  10  psi,  V  =  0.3 

Finite  Element  Breakdown:  20  elements  in  the  first  2  inches 

from  the  frse  end,  15  alements  in  the  next  3  inches j  and 
10  elements  in  the  last  5  inches 

Ring  Leading  Q  at  s 
(lb/in)  (in) 

1.0  0 

0.75  0.3 

0.50  0.6 

0.25  0.9 

(b)  Sandwich  Conical  Shell 

Same  as  in  (a) ,  except 

Thickness:  0.01  in.  for  top  and  bottom  layer  and  0.08  in.  for 
the  core 

7 

Material:  E  =  10  psi,  V  =  0.3,  isotropic  for  cover  layers,  and 
4  6 

E  =  10  psi,  'J  =  0. 3 ,  g  =  10  psi  for  the  core 

(c)  Two-Layer  Conical  Shell 

Same  as  in  (a) ,  except 

Thickness:  0.5  in.  for  cover  material  and  0.5  in.  for  base  material 

6 

Material:  Glass  Phenolic  cover,  L  =  10  psi,  v  =  0.3,  isotropic 
Aluminum  base,  E  =  10^  psi,  V  =  0.3  isotropic 
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TABLE  « 


PERCENTAGE  ERROR  OF  THE  LOWEST  FREQUENCY  PREDICTED  FOR 
A  SIMPLY-SUPPORTED  SQUARE  THIN  PLATE  USING  ELEM1 

W  -  W 

( - -  -  1)  x  100 

(i) 

o 


MESH+ 

DOF 

LUMPED  MASS 

HR  MASS 

2x2 

27 

-  5.2 

5.1 

4x4 

75 

-  1.3 

1.2 

6x6 

147 

-  0.6 

0.6 

8x8 

243 

-0.4 

0.3 

2  4 

Thickness  =0.1  in..  Side  length  -  20  in.,  E  =  1  psi,  V  =  0.3,  P  =  1( lb-sec  )/in 
»  Exact  Lowest  Frequency  (Ref.  40)  =  0.00149  rad/sec 

+Uniform  mesh  in  a  quarter  of  the  plate 

TABLE  9 

i 

PERCENTAGE  ERROR  OF  THE  LOWEST  FREQUENCY  PREDICTED  FOR 
A  CLAMPED  RECTANGULAR  THIN  PLATE  USING  ELEM1 
W  -  OJ 

( - -  -  1)  x  100 

u 


+ 

MESH 

DOF 

LUMPED  MASS 

HR  MASS 

2x2 

27 

-  7.1 

4.0 

4x4 

75 

-  2.0 

1.1 

6x6 

147 

-  1.0 

0.3 

r>..o 

A  “I 

OAV 

i 

\J  •  ‘  » 

Thickness  =0.1  in..  Short-side  length  =  20  in..  Long-side  length  =  40  in., 
E  =  1  psi,  V  =  0.3;  p  =  1 (lb- sec2) /in4 

tj  =  Exact  Lowest  Frequency  (Ref.  40)  =  0.00187  rad/sec 
o 

+Umform  mesh  in  a  quarter  of  the  plate. 
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TABLE  10 


FREQUENCY  PREDICTIONS  FOR  THE  FIRST  SEVEN  SYMMETRIC  MODES  OF  VIBRATION 
OF  A  SIMPLY-SUPPORTED  RECTANGULAR  THIN  PLATE  USING  ELEM1  AND  THE 

HYBRID- RATIONAL  HASS  MATRIX 


MODE* 

m-n 

EXACT** 

U) 

6x12 

MESH* 

3x8  MESH++ 

z 

Error  (%) 

u 

ERROR  (j) 

1-1 

1.25 

1.2565 

0.52 

1.2540 

0.32 

3-1 

3.25 

3.2363 

1.12 

3.3340 

2.58 

5-1 

7.25 

7. 5021 

3."74 

7.9220 

9.27 

1-3 

9.25 

9.7807 

5.74 

9.5198 

2.92 

3-3 

11.25 

11.7754 

4.40 

11.5685 

2.83 

7-1 

13.25 

14.2323 

7.41 

16.0343 

21.01 

5-3 

15.25 

15.9259 

4.43 

_ 

16.545 

5.28 

Mode  Shape:  sinCnTr/a'  sin  {nn/b} 

**  2  ) -  2  2  2  2 

Ref.  40  exact  solution;  w  =  tt  /D/?h  (m  /a  +  n  /b  )  ;  a  =  2b; 

normalized  (i)  =  {m^/4  +  n^)  ;  D  -  Eh^/[12(1  -  V^)]. 

Results  of  a  6x2  uniform  mesh  in  a  quarter  of  the  plate  (6  in  the  long 
dimension);  total  DOF  =  273;  total  computing  time  on  IBM  370/155  =  5.9  min. 

++Results  of  an  8x8  uniform  mesh  in  a  quarter  of  the  plate;  total  DOF  =  243; 
total  computing  time  on  IBM  370/155  =  6.6  min. 
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TABLE  11 


PERCENTAGE  ERROR  OF  THE  LOWEST  FREQUENCY  PREDICTED  FOR  A 
TWO- LAYER  CROSS-PLY  (0°/90°)  SQUARE  THIN  PLATE  USING 
ELEKR  AND  THE  HR  MASS  MATRIX 


10  - 

CO 

o 

1)  x  100 

(  <0 

p 

MESH 

DOF 

HR  MASS 

2x2 

45 

3.0 

4x4 

125 

-  0.5 

6x6 

245 

-  0.13 

8x8 

405 

-  0.16 

VE2  ’  G12/E2  ■  °'5'  V12  *  °'25'  C  '  1 

(jJq  =  Lowest  frequency  calculated  in  Ref.  22  =  0.00614  rad/sec 


TABLE  12 

DATA  FOR  THE  VIBRATION  ANALYSIS  OF  A  SIMPLY -SUPPORTED  TWO-LAYER 
CROSS-PLY  (0°/90°>  INFINITE  STRIP 


(a)  Material  Properties 

Bottom-Layer  =  E.  =  1.0922,  E  =  41.983,  V  =  0.00520,  v  =  0.2,  G  =  G 
L  &  a  j 

Top-Layer  =  Ex  =  43.516,  E2  =  1.0455,  v12  =  0.20812,  V23  =  0.2,  G12  =  =  1 


(b)  Finite-Element  Discretization 

Eight  ELEMZ  elements  in  half  strip  (£/2) and  two  sublayers  for  each  layer; 
i.e.,  treating  the  two- layer  plate  as  a  four-layer  plate. 
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TABLE  13 


COMPARISON  OF  PREDICTIONS  FOE  THE  LOWEST  FREQUENCIES  OF  3-LAYER 
SIMPLY-SUPPORTED  SQUAPE  PLATES 


CASE 

VS 

FREQUENCY (A) 

[  ERROR  % 

CPT 

ELEMZ+ 

EXACT 

CPT 

ELF.MZ 

1 

1 

1 

0.077054 

0.0745 

3.4 

1.1 

2 

1 

2 

0.093994 

0.090973 

0.089986 

4,5 

i.l 

3 

1 

5 

0.13239 

0.12437 

0.123072 

7-JU 

1.1 

4 

1 

15 

0.215642 

0.18529 

0.103664 

17.4 

0.9 

5 

2 

15 

0.196852 

0.16897 

0.167574 

17.5 

0.8 

6 

3 

15 

0.18225 

0.15633 

0.155082 

17.5 

0.8 

Isotropic,  =  V2  =  V3  =  0-3'  =  G3'  ^  =  Mass  density,  d^  =  d^ 

h^/a  =  0.01,  h^/a  = 

A  =  w (d^h^/G^) A/2 (  w  =  angular  frequency  (rad/sec) 

CPT  =  classical  plate  theory,  Ref.  42 
EXACT  =  Ref.  42 

Finite-Element  Discretization  -  uniform  4x4  mesh  in  a  quarter  of  the  plate 
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TABLE  14 

PROPERTIES  OF  A  THREE-LAYER  CYLINDRICAL  SHELL 


LAYER 

THICKNESS 
(in. ) 

C  * 
11 

C22 

C33 

C12 

C13 

C23 

C44 

C55 

C66 

Inside 

0.2 

33.03 

3.032 

33.32 

0.3998 

10.03 

1.03 

1.154 

1.154 

1.154 

Middle 

0.5 

33.02 

3.032 

RK9 

0.3998 

10.03 

1.03 

11.54 

0.1154 

0.1154 

Outside 

0.3 

33.02 

3.032 

0.3998 

10.03 

1.03 

1.154 

1.154 

1.154 

*  6 

Elastic  constants  C.  .  defined  by  0.  =  C..  E.,  i.  j  =  1,  2  ,...6  (units  in  10  psi) 

n  i  13  3  '  ' 


TAELE  15 

COMPARISON  OF  PREDICTIONS  FOR  THE  LOWEST  FREQUENCIES  OF  1 -LAYER 
SIMPLY-SUPPORTED  CYLINDRICAL  SHELLS  IN  AXISYMMETRIC  VIBRATION 


ti/Z 

a 

(in) 

FREQUENCY (RAD/SEC) 

ERROR* 

%) 

CST 

RST 

ELEMR"*' 

ELEMZ+ 

ET 

CST 

RST 

ELEMR 

ELEMZ 

0.01 

0.5 

171.65 

171.65 

■Bli 

172.91 

171.65 

0.0 

0.0 

0.7 

0.7 

0.5 

177.54 

177.28 

■sa 

177.35 

177.35 

0.1 

0.0 

0.0 

0.0 

0.1 

233.14 

228.15 

230.10 

228.63 

2.0 

-  0.2 

0.6 

0.2 

0.20 

.05 

647.82 

560.22 

593.67 

574.00 

568.52 

13.7 

-  1.6 

4.4 

0.9 

0.40 

.025 

2502.9 

1607.8 

1074.2 

1728.5 

1704.4 

46.8 

-  5.7 

9.9 

SB 

0.60 

.016 

5620.8 

2751.2 

_ 

3411.4 

3099.1 

3052.3 

84.1 

-  9.9 

11.7 

■B 

.012 

9989. 3 

3886.9 

5003.0 

4557.3 

4495.2 

122.2 

-13. 5 

11.3 

1.4 

39 

.01 

L5607.0 

5005.5 

6593.2 

6062.4 

5991.7 

160.5 

-16.5 

10.0 

1.2 

H/£  ^Thickness- to- length  ratio.  Mean  Radius  R  =  10 
a  =  Width  of  flat-plate  element 

CST  =  Classical  Shell  Theory,  Ref.  50. 

rpv,„„„.,  Onf  ,10 

IWi  —  J  ,  *M-.  *  ..  * 

ET  =  Elasticity  Theory,  Ref.  48 
+The  HR  mass  matrix  is  employed 
++ERR0R  =  Compared  with  the  Elasticity  Solution  of 


Ref.  48. 
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TABLE  16 


DATA  FOR  A  FREE  SINGLE-LAYER  ISOTROPIC  THIN  RECTANGULAR 
PLATE  SUBJECTED  TO  THERMAL  LOADING 

(a)  Material  Properties  : 

E  =  10.4  x  10°  psi ,  V  =  0.3,  a  -  12.7  x  10  in/in-°F 

(b)  Dimensions: 

Rectangular  plate  with  long-side  length  =  2a  =  36  in,  short-side 
length  =  2b  =  24  in,  and  thickness  =  0.25  in 

(c)  Temperature  Distribution 

Linear  variation  in  the  short-span  direction  and  constant 
in  the  long-span  direction.  Maximum  temperature  rise  T^ 
is  150  °F  (see  Fig.  28b) 

(d)  Finite  Element  Discretization: 

Uniform  5x4  mesh  in  a  quarter  of  the  plate  with  ELEMZ 
"degeneralized”  for  this  single-layer  problem 

TABLE  17 

DATA  FOR  A  CLAMPER  RECTANGULAR  SINGLE -LAYER  ISOTROPIC  THIN  PLATE 
SUBJECTED  TO  THERMAL  LOADING 

(a)  Material  Properties: 

E  =  1.0  psi,  V  =  0.3,  a  =  1.0  x  10  6  in/in-°F 

(b)  Dimensions: 

Rectangular  plate  with  long-side  length  =  32  in,  short-side  length 
=  16  in  and  thickness  =  0.16  in 

(c)  Temperacure  Distribution: 

Parabolic  in  the  short-span  direction  and  constant  in  the 
long-span  direction  (see  Fig.  .id) 

(d)  Finite-Element  Discretization: 

Uniform  8x0  mesh  in  a  quarter  of  the  plate  with  ELEMZ 
"degeneralized"  for  tiiis  single-layer  problem 
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TABLE  Id 


DATA  FOR  A  SIMPLY-SUPPORTED  TWO-LAYER  CROSS-PLY  (0°/90°) 
INFINITELY-LONG  THIN  STRIP 

(a)  Material  Properties: 

Vb2  ‘  25'  °12/E2  *  °-5'  G23/E2  ‘  °‘2'  "  S3  "  »•“* 

E  =  1  psi  (Longitudinal  a) /Transverse  a)  =  1/3, 

*  -6 
Longirudinal  a  =  1.0  x  10  in/in-°F 

(b)  Dimensions: 

Short-span  length  =  16  in.  Thickness  =  0.16  in  (0.08  in.  for 
each  layer) 

(c)  Temperature  Distribution: 

1  6F  uniformly  distributed  over  the  entire  strip 

(d)  Finite-Element  Discretisation: 

Eight  ELEMZ  elements  in  a  half  span 

(e)  Boundary  Condition:  only  transverse  displacements  along 

the  supported  edges (x  =  +  8  in.)  are  constrained. 
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TABLE  19 

DATA  FOR  A  SIMPLY -SUPPORTED  THREE-LAYER  CROSS-PLY  (0°/90o/0!>) 
CYLINDRICAL  SHELL  SUBJECTED  TO  THERMAL  LOADING 

(a)  Material  Properties: 

E/E  =  25,  C  /E  »  0. 5;G_-/E  =  0.2,  V  =  V  =  0.25, 

12  12  2  23  2  12 

E  =  1  psi  (Longitudinal  a) /(Transverse  a)  =  1/3,  Longi- 
tudinal  a  =  1.0  x  10  in/in-°F 

(b)  Dimensions: 

Mean  Radius  =  20  in.,  Thickness  =  1  in.  (0.2  in.  for  inside 
layer,  0.5  in.  for  middle  layer  and  0.3  in.  for  outside 
layer).  Total  length  =  16  in. 

(c)  Temperature  Distribution: 

1  °F  uniformly  distributed  over  the  entire  shell 

(d)  Finite-Element  Discretization: 

Eight  ELEMZ  elements  in  half  strip 

(e)  Boundary  Conditions:  only  radial  displacements  along  the 

supported  edges  are  constrained. 
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ERROR*  IN  v  (PER  CENT) 


0  100  200  300  400  500  600 


QUARTER  PLATE  DOF 

FIG.  2  CENTRAL  DEFLECTION  OF  A  SIMPLY -SUPPORTED,  CENTRALLY- 
LOADED,  SQUARE,  ISOTROPIC,  THIN,  SINGLE-LAYER  PLATE 
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ERROR*  IN  W  (PER  CENT) 


PIG.  3  -CENTRAL  DEFLECTION  OF  A  SIMPLY -SUPPORTED,  UNIFORMLY- LOADED, 
SQUARE,  ISOTROPIC,  THICK,  SINGLE-LAYER  PLATE 
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NONDIMENSIQNAL  CENTRAL  DISPLACEMENT.  W  D/Pa 


PRESENT  (ELEM1) 


EQUILIBRIUM  MODEL  (REF.  19) 


HYBRID  STRESS  MODEL.  STRESS  BOUNDARY 
CONDITIONS  ENFORCED  (REF.  16) 


-  ANALYTICAL  (REF.  21) 


□  DISPLACEMENT  MODEL  (REF.  20) 


•  HYBRID  STRESS  MODEL,  STRESS 
BOUNDARY  CONDITIONS  NOT 
ENFORCED 


100  200  300  400  500  600  700  800  900 


LX iF  FOR  ENTIRE  PLATE 


FIG.  7  CENTRAL  DEFLECTION  OF  A  SKEWED,  SINGLE-LAYER,  ISOTROPIC 
THIN  PLATE  UNDER  UNIFORM  LOAD 


QUARTER  PLATE  DOF 

FIG.  8  CENTRAL  DEFLECTION  OF  A  CLAMPED,  TWO-LAYER  RECTANGULAR, 
CROSS-PLY  PLATE  UNDER  UNIFORM  LOAD 


0  100  200  300  400  500  600  700  800  900  1000 

DOF  FOR  ENTIRE  PLATE 


FIG.  9  CENTRAL  DEFLECTION  OF  A  CLAMPED,  TWO-LAYER,  SQUARE,  CROSS-PLY 
PLATE  UNDER  UNIFORM  LOAD 
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50 


£  30 


*  20 
C« 


ELEMR 


COMPARED  WITH  CPT 
SOLUTION  (REF.  23) 


(+  45*)  and  (0*/90°) 


(+  45*) 


(0*/90# ) 


8 


DOF  FOR  ENTIRE  PLATE 

FIG.  10  CENTRAL  DEFLECTION  OF  A  CLAMPED,  TWO- LAYER,  SQUARE 
ANGLE-PLY  PLATE  UNDER  UNIFORM  LOAD 


ERROR*  IN  %»  (PER  CENT) 


nmwwaw!  wppw 


LEGEND  FOR  FIG. 12 

-  ELASTICITY  (REF.  1) 

-  PRESENT  (ELEMZ) 

- FINITE  ELEMENT  (REF. 5) 

-  CPT  (REF.  1) 

(SEE  TABLE  4a  FOR  DIMENSIONS  ANU 
NONDIM.  (~)  QUANTITIES) 


•COMPARED  WITH  EXACT  SOLUTIONS  (REF.  1) 

(a)  Central  Displacement,  (b)  In-plane  Displacement,  u 


(c)  Normal  Stress,  Ox 


(d)  Shear  Stress,  0 

xz 


FIG.  12  DISPLACEMENT  AND  STRESS  SOLUTIONS  FOR  AN  INFINITE  THREE-LAYER 
CROSS-PLY  PLATE  WITH  SIMPLY -SUPPORTED  EDGES,  SUBJECTED  TO 
SINUSOIDAL  LOADING 
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LEGEND  FOR  FIG. 13 


-  ELASTICITY  (REF.  2) 

-  CPT  (REF.  2) 

-  PRESENT  (ELEMZ) 

0  FINITE  ELEMENT  (REF.  6) 

(See  b2  of  Table  4  for  Dimensional  and  Non- 
Dimensional  (  )  Quantities) 


*/h 


*/h 


(c)  Transverse  Shear,  a 

xz 


z/h 


FIG.  13  DISPLACEMENT  AND  STRESS  SOLUTIONS  FOR  A  SIMPLY-SUPPORTED 
THREE-LAYER  CROSS-PLY  SQUARE  PLATE  SUBJECTED  TO  SINU¬ 
SOIDAL  LOADING 
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CENTRAL  RADIAL  DEFLECTION  V  X  10  (IN) 


100 


i  r 


-  SABOR  4  (35  ELEMENTS,  REF.  28) 

-  PRESENT  (35  ELEMENTS,  ELEMl) 

°  PRESENT  (45  ELEMENTS,  ELEMl) 

THICKNESS  *  0.025  IN 


_ I _ I _ I _ L 

2  3  4  5 

MERIDIONAL  DISTANCE  F RDM  THE  FREE  END  (IN) 

(b)  Radiax  Displacement  w 

FIG.  17  CONCLUDED 
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MERIDIONAL  MOMENT  x  100  (IN-LB) 


0  1  2  3  4  5 


MERIDIONAL  DISTANCE  FROM  THE  FREE  END  (IN) 

(a)  Meridional  Moment 

FIG.  18  MERIDIONAL  MOMENT  AND  RADIAL  DISPLACEMENT  OF  A  THIN,  ISOTROPIC 
SINGLE-LAYER,  CONICAL  SHELL  UNDER  RING  LOADS 
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0 


1 


4 


5 


10 


2  3 

MERIDIONAL  DISTANCE  PROM  THE  FREE  END 
(b)  Radial  Displacement,  w 

FIG.  18  CONCLUDED 
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MERIDIONAL  DISTANCE  FROM  THE  FREE  END  (IN) 

FIG.  15  MERIDIONAL  MOMENT  OF  TWO  MODERATELY  THICK,  SINGLE-LAYER, 
ISOTROPIC',  CONICAL  SHELLS  CINDER  RING  LOADS  ‘ 


RADI  Ail  DISPLACEMENT,  W  X  10  (IN) 


% 

•\ 

I  — 

\ 

-  \ 

\ 

\ 

■  \  - 
I 

\ 

\  \ 

\  \ 

■  \  \ 

\  > 

.  w 

\  » 

\  \\ 

'  \  ''  \ 
\  \  \\ 


-  THICKNESS  -  0.5  IN.  PRESENT,  ELEMl 

(SINGLE-LAYER) 

-  THICKNESS  -  1.0  IN,  PRESENT,  ELEMl 

(SINGLE-LAYER) 

-  THICKNESS  -  1.0  IN,  SABOR  4  (REF. 28) 

(SINGLE-LAYER) 


THICKNESS  *  1.0  IN,  PRESENT,  ELEMR 
(TWO-LAYER) 


\  \\\ 
\  \  '  > 


MERIDIONAL  DISTANCE  FROM  THE  FREE  END  (IN) 

FIG.  20  RADIAL  DISPLACEMENTS  OF  TWO  SINGLE-LAYER  AND  ONE  TWO-LAYER 
CONICAL  SHELLS  UNDER  RING  LOADS 


10  (IN) 


meridional  distance  from  the  free  END  (IN) 

(a)  Maximum  Meridional  Normal  Straaaea  in  tha  Cover  and  in  the  Bare 

FIG.  22  MAXIMUM  MERIDIONAL  AND  CIRCUMFERENTIAL  NORMAL  STRESSES  IN  THE  COVER  AND 
IN  THE  BASE  OF  A  TWO- 1 X YE R  CONICAL  SHELL  UNDER  RING  LOADS 


(ISJ) 


'S3SS3HXS  TWHON  lVUN3H3iWnCMID  WnWIXVW 
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ERROR*  IN  THE  LOWEST  FREQUENCY,  u  (PER  CENT) 


FIG.  23  COORDINATE  TRANSFORMATION  FOR  A  QUADRILATERAL  ELEMENT 


FIG.  24  THE  LOWEST  FREQUENCY  OF  A  SIMPLY-SUPPORTED,  SINGLE-LAYER 
THIN,  ISOTROPIC,  SQUARE  PLATE 
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FIG.  25  THE  LOWEST  FREQUENCY  OF  A  CLAMPED,  SINGLE-LAYER,  ISOTROPIC, 
THIN,  RECTANGULAR  PLATE 
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T  FREQUENCY,  U)L(0/10) 


C6 


U 


FIG.  26  THE  LOWEST  FREQUENCY  OF  A  CLAMPED,  TWO-LAYER, 
CROSS-PLY  THIN,  SQUARE  PLATE 


10 


8 


7 

6 

5 

4 


a 

a  2 


0 


0  0.4  0.8  1.2  1.6  2.0 


THICKNESS-TO-SPAN  RATIO,  h/L 

FIG.  27  THE  LOWEST  FREQUENCIES  OF  INFINITE  TWO-LAYER,  CROSS-PLY 
PLATES  WITH  SIMPLY -SUPPORTED  EDGES 
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LINEAR  IN  X  AND  y 


7/ 


f/7 


1^2  L_ 


AT  (x,y,z) 


Ith  LAYER 


R  /  PARABOLIC  IN  Z 


(a)  Corner  Points  and  Temperature  Assumptions 


T  +  xT^  +  yT. 


*4.  +-  ■*  t 


TOP  SURFACE 


-T*- 


T  +  xT  +  yT 
^  ml  m2  m3 

- X*  middle  surface 


.T.  +  xT  +  yT 

Y  bi  2  3 


BOTTOM  SURFACE 


(b)  Temperature  Assumptions  on  Three  Surfaces 

FIG.  28  TEMPERATURE  ASSUMPTIONS  FOR  ONE  LAYER  OF  A  MULTILAYER 
QUADRILATERAL  ELEMENT 
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(a)  Dimensions  and  Coordinate  System 


FIG.  29  DIMENSION,  COORDINATE  SYSTEM,  AND  TEMPERATURE  DISTRIBUTION 
OF  A  FREE,  SINGLE-LAYER,  THIN,  ISOTROPIC  RECTANGULAR  PLATE 


ORIGINAL  SHAPE 
DEFORMED  SHAPE 


Q.QljIN.  SCAU;  FQR  displacement 


FIG.  30  DISPLACEMENTS  OF  A  QUARTER  OF  THE  PLATE  OF  FIG.  29 
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t  »  7.2  IN 

4  8  12 


DISTANCE,  y  (IN) 


(b) 


Transverse  Normal  Stresses ,  0 

y 


•  THEORY 
(REF.  53) 


°  EXPERIMENT 
(REF.  53) 


•  PRESENT 
(ELEMZ) 


FIG.  31  CONTINUED 
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THEORY 

0  4  B  12 

DISTANCE,  y  (IN) 

(c)  Shear  Stresses,  O 

xy 

FIG.  31  CONCLUDED 


THEORY 
(REF.  53) 


O  EXPERIMENT 
(REF.  53) 

•  PRESENT 
(ELEMZ) 


r 


DEFLECTION 


8 


10  (IN) 


(a)  pinna iona  and  Coordinate  System 


012345678 


DISTANCE,  x  (IN) 

(b)  Deflection,  w 


10  PSI 

! — I 


x-.5  x-1.5  x-2.5  x«3. 5  x-4.5  x*5.5  x-6.5  x«7.5  (IN) 

NORMAL  STRESSES,  0x  (PSI) 

(c)  Normal  Streases,  0 

FIG.  33  DIMENSIONS,  COORDINATE  SYSTEM,  DEFLECTIONS,  AND  STRESSES  OF 
AN  INFINITE  TWO-LAYER  CROSS-PLY  THIN  PLATE  SIMPLY  SUPPORTED 
ALONG  EDGES  AND  SUBJECTED  TO  UNIFORM  TEMPERATURE 


m 


IPWIIIWPSP 


(a)  Radial  Deflections,  w 


(D)  Axial  Deflections,  u  (in) 


PIC.  34  RADIAL  AND  AXIAL  DEFLECTION  OP  A  SIMPLY-SUPPORTED  THREE-LAY^ 
CROSS-PLY  CYLINDRICAL  SHELL  SUBJECTED  TO  IwSSS' T^T^f 
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APPENDIX 


COMPUTATIONAL  SCHEMES  FOR  STATIC  AND  VIBRATION  ANALYSIS 


Since  the  purpose  of  this  study  is  to  develop  and  evaluate  new  useful 
finite  elements,  example  problems  chosen  are  often  relatively  simple  in  geo¬ 
metry  so  that  results  can  be  compared  with  those  of  other  existing  solutions 
and  computation  can  be  completed  economically.  Consequently,  the  computa¬ 
tional  schemes  chosen  for  use  in  this  study  are  suitable  for  medium-sized  prob¬ 
lems  with  total  degrees-of-freedom  no  more  than  about  1000. 

A  thermal  stress  analysis,  as  shown  in  Section  4,  can  be  reduced  to  a 
static  analysis  by  transforming  the  initial  thermal  strains  into  an  equivalent 
loadii,.;  vector.  Thus,  for  both  static  and  thermal  stress  analyses,  the  system 
equations  to  be  solved  may  be  represented  by 


k  r=F 


(A.  1) 


where  £  is  the  system  stiffness  matrix,  q*  is  the  unknown  nodal  displacement 
vector,  and  F  is  the  force  vector.  To  solve  Eq.  A.l  the  method  of  triangular 
factorization  is  used.  The  £  matrix  is  decomposed  into  three  matrices: 


K=ULr 


(A.  2) 


where  L  is  a  lower  triangular  matrix  with  zero  elements  in  the  upper  triangle 
and  unity  on  the  diagonal,  and  D  is  a  diagonal  matrix.  Note  that  the  decom¬ 
position  is  easily  achieved  by  successive  substitution  starting  from  the  first 
row-and- column  element  of  JC.  Combining  Eq.  A.l  and  Eq.  A. 2,  one  obtains 


L&L^r-F 


Defining  a  new  unknown  vector 


A.  'V  V  -A/ 

one  obtains 


(A. 3) 


(A. 4) 


L  Y=F 


(A.  5) 
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Now  Eq.  A. 5  can  be  solved  readily  for  JC  since  l  is  a  lower  triangular  matrix; 
the  solution  of  X  can  be  obtained  by  successive  substitution  starting  from  the 
first  unknown.  Once  V  is  obtained,  q*  can  be  obtained  similarly  from  Eq.  A. 4 

T  ~ 

by  noting  that  D  L  is  U)  upper  triangular  matrix.  Thus,  two  steps  a..- 
volved  in  the  solution  of  Eq.  A.l:  (a)  decomposition,  which  also  involves  only 
successive  substitution,  and  (b)  successive  substitution  to  obtain  X  and  then 

For  vibration  analysis,  a  system  mass  matrix  is  added  and  the  system 
equations  are,  from  Eq.  3.15, 

(K  (A. 6) 

where  u  is  the  natural  frequency  to  be  determined. 

Many  eigenvalue  schemes  can  be  used  to  solve  Eq.  A. 6,  but  it  is  con¬ 
venient  in  programming  to  make  use  of  the  triangular  factorization  scheme  of 
the  static  analysis.  Thus,  the  eigenvalue  solution  using  the  Sturm-sequence 
method  [42,43,44]  is  selected  for  this  study.  The  Sturm  sequence  property, 
when  applied  to  Eq.  A. 6,  states  that  for  any  given  value  of  uj,  the  number  of 

successive  sign  changes  in  the  determinants  of  the  leading  minors  of  matrix 
2 

[K  -  uj  £}J  equals  the  number  of  eigenvalues  that  are  less  than  w.  Obviously, 
this  property  can  be  utilized  to  locate  or  "bracket"  any  eigenvalue  by  simply 
assuming  a  series  of  oi's  and  counting  the  number  of  sign  changes  in  the  de¬ 
terminants. 

An  important  observation  from  the  decomposition 

-  L  I>LT  (A. 7) 

is  that  the  number  of  sign  changes  of  the  determinant  of  the  leading  minors 
equals  the  number  of  negative  values  in  the  diagonal  of  D.  Therefore,  it  is  a 
simple  task  to  determine  the  number  of  eigenvalues  less  than  a  given  w  once  the 
decomposition  of  Eq.  A. 7  is  complete. 

Using  the  above  procedure,  any  eigenvalue  can  be  isolated,  and  lower  and 
upper  bounds  can  be  established.  But  for  accurate  determination, it  is  necessary 
to  use  an  interpolation  scheme  to  reduce  further  the  bracket  interval  of  an 
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eigenvalue  to  desirable  accuracy.  In  this  study,  a  parabolic  curve  of  co  vs. 

2 

the  determinant  of  ( K  -  w  M)  is  fitted  between  the  bounds  of  an  eigenvalue, 
and  the  interpolated  value  of  td  is  the  estimated  eigenvalue.  This  interpola¬ 
tion  process  is  continued  until  changes  in  id  is  within  a  certain  prescribed 
amount  of  tolerance  (for  example,  1/100  of  the  calculated  w) . 

Once  an  eigenvalue  is  obtained,  the  corresponding  eigenvector  is  ob¬ 
tained  by  solving  Eq.  A. 6  with  one  of  the  degrees-of-freedom  constrained. 
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